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ABSTRACT 


The  engineering  problem  is  the  design  of  submerged  minimum  drag 
axlsymmetric  vehicles  for  a  specified  enclosed  volume  and  constant 
speed.  Propulsion  is  not  considered  so  that  drag  rather  than  power  Is 
to  be  minimized.  Drag  reduction  Is  to  be  accomplished  solely  through 
manipulation  of  the  vehicle  shape;  other  means  of  drag  reduction,  such 
as  polymer  injection  into  the  boundary  layer,  are  not  considered. 

The  optimization  problem  is  formulated  as  a  ncngradlent  search  in 
a  finite  constrained  parameter  space.  Two  classes  of  bodies,  described 
by  five  and  eight  parameters,  are  considered.  The  bodies  are  con¬ 
strained  to  be  well  behaved  based  on  previous  hydrodynamic  experience. 
The  drag  model,  valid  for  nonseparating  flows,  consists  of  computer 
programs  available  In  the  literature  and  is  representative  of  state-of- 
the-art  drag  prediction  methods.  The  requirement  for  nonseparating  flow 
represents  an  additional  constraint  on  the  optimization  problem.  Two 
optimization  methods  representing  diverse  search  philosophies  are  used 
to  obtain  the  optimal  solutions.  These  Include  Box's  Complex  Method 
and  Powell's  Method  of  Conjugate  Directions  used  In  conjunction  with 
a  penalty  function. 

The  results  show  that  significant  drag  reduction  is  possible 
through  shape  manipulation.  Reductions  of  one-quarter  to  one-third 
below  the  best  existing  designs  have  been  obtained.  All  optimal  designs 
exploit  laminar  boundary  layers.  If  laminar  flow  Is  not  allowed,  then 
drag  reduction  below  the  best  existing  designs  apparently  must  be 
accomplished  by  means  other  than  shape  manipulation.  The  optimal  lami¬ 
nar  shape  is  a  strong  function  of  Reynolds  number,  ranging  from  quite 


slender  at  low  values  to  quite  "fat”  at  high  values.  The  minimum  drag 
shapes  have  a  high  sensitivity  to  early  transition.  Suboptlmal  low 
drag  bodies  without  this  characteristic  are  used  for  hydrodynamic  de¬ 
sign.  Extensive  runs  at  one  Reynolds  number  suggest  that  there  Is  a 
unique  global  minimum  drag  shape  at  each  Reynolds  number. 
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CHAPTER  1 
INTRODUCTION 


Of  fundamental  interest  in  the  field  of  fluid  mechanics  is  the  study 
of  fluid  forces  exerted  on  a  moving  body.  In  the  area  of  hydrodynamics, 
a  time-honored  goal  of  the  naval  architect  has  been  the  reduction  of 
vehicle  resistance  by  means  of  vehicle  shape  as  well  as  by  more  elaborate 
schemes  [1],  such  as  viscoelastic  polymer  injection  into  the  boundary 
layer. 

We  are  considering  here  the  shape  and  resulting  resistance,  or  drag, 
of  bodies  submerged  in  an  incompressible  fluid.  The  applications  are 
directed  toward  vehicles  such  as  torpedoes  and  submarines  deeply  sub¬ 
merged  in  water.  The  primary  function  of  such  vehicles  is  to  provide  an 
enclosed  volume  in  which  a  payload  is  carried. 

The  purpose  of  the  present  study  is  to  include  the  hydrodynamics 
as  part  of  an  optimization  problem  in  which  the  vehicle  drag  is  to  be 
minimized  through  shape  manipulation.  Other  means  of  drag  reduction  are 
not  considered.  The  hydrodynamic  problem  is  made  more  tractable  by  re¬ 
stricting  the  analysis  to  the  class  of  axisymmetric  vehicles  (bodies  of 
revolution)  without  appendages  immersed  in  ax i symmetric  flow  (zero  angle 
of  attack). 

Stated  somewhat  more  precisely,  the  optimization  problem  involves 
the  -following:  for  a  specified  incompressible  fluid,  vehicle  volume, 
an i  constant  vehicle  speed,  find  the  axisymmetric  body  shape  which 
minimizes  the  drag.  The  vehicle  whose  shape  minimizes  the  drag  is 
defined  as  the  optimum  body.  The  equivalent  nondimensional  optimization 
problem  may  be  stated  as  follows:  for  a  specified  Reynolds  number 

i 

Ry  -  V3/\>,  find  the  vehicle  shape  which  minimizes  the  drag 
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coefficient  Cg  =  D/(j.p  U£  VT),  where  U»  is  the  constant  vehicle 
speed,  V  is  the  vehicle's  enclosed  volume,  v  is  the  fluid  kinematic 
viscosity,  D  is  the  vehicle  drag,  and  p  is  the  fluid  density. 

The  primary  concern  here  is  the  design  of  the  largest  vehicle  volume 
with  the  lowest  possible  drag  at  a  given  speed;  there  are  no  prescribed 
constraints  on  the  body.  From  the  hydrodynamics  analysis  point  of  view, 
the  actual  volume  and  velocity  separately  are  not  important;  the  body 
shape  and  characteristic  Reynolds  number  Rv  completely  establish  the 
fluid  dynamics.  Thus  a  particular  body  shape  has  the  same  drag  coeffi¬ 
cient  Cq  over  a  wide  range  of  velocities  and  volumes  so  long  as  the 
Reynolds  number  Rv  is  unchanged. 

Other  design  problems  may  be  more  conveniently  handled  by  using  a 

i 

characteristic  Reynolds  number  based  on  something  other  than  (volume)7. 

For  example,  the  torpedo  must  have  a  fixed  maximum  diameter  (constant 
frontal -area);  it  is  more  convenient  to  specify  a  Reynolds  number  based 
on  maximum  diameter  rather  than  volume.  Other  characteristic  lengths 

which  may  be  useful  in  design  problems  include  body  length  and  (wetted 

1 

area)2.  It  is  emphasized  that  the  optimization  formulation  in  the  pre¬ 
sent  study  applies  to  any  of  these  design  problems. 

It  is  not  known  if  a  unique  optimum  body  exists  for  each  Reynolds 
number.  That  is,  there  may  be  an  entire  class  of  bodies  which  have  the 
same  minimum  drag  coefficient  at  a  given  Reynolds  number.  If  such  is 
the  case,  then  the  designer  must  introduce  additional  considerations  or 
constraints  to  obtain  the  one  design  best  suited  for  his  application. 

For  example,  the  optimization  problem  ignores  the  sensitivity  of  drag 
changes  with  body  velocity  variations  and  ignores  the  effect  of  angle 
of  attack  on  drag.  If  the  designer  has  several  minimum  drag  body  de¬ 
signs,  he  considers  these  kinds  of  ideas  in  making  the  final  selection. 
Such  considerations  are  outside  the  scope  of  the  present  study.  A  more 
complete  discussion  on  design  constraints  is  given  in  Chapter  3. 

The  application  of  formal  optimization  methods  to  the  drag  minimiza¬ 
tion  problem  has  not  appeared  in  the  literature  to  date.  One  reason  for 
this  absence  is  the  fact  that  no  reliable  drag  model  for  arbitrary  axi- 
symmetric  bodies  in  incompressible  flow  is  available.  More  fundamentally. 
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there  is  an  incomplete  understanding  of  the  fluid  flow  physics  in  the 
boundary  layer  transition  region,  in  the  turbulent  boundary  layer  develop¬ 
ment,  in  the  turbulent  boundary  layer  separation  region,  and  in  the  wake 
region  following  separation.  All  of  these  phenomena  are  of  primary  Im¬ 
portance  in  drag  prediction.  In  addition,  there  are  other  factors  which 
complicate  the  prediction  of  drag  in  practice  [2],  such  as  the  effects 
of  ambient  turbulence,  body  surface  waviness,  and  body  vibration  on  the 
boundary  layer  development. 

With  the  lack  of  reliable  drag  models,  all  published  work  to  date 
directed  toward  the  design  of  low  drag  axi symmetric  bodies  in  incompress¬ 
ible  flow  has  been  mainly  experimental  in  nature.  Two  studies  will  be 
noted  here. 


Series  58  Stud> 


An  experimental  study  of  drag  for  a  systematic  series  of  axisym- 
metric  bodies  was  reported  by  Gertler  in  1950  [3].  The  series,  known 
as  "Series  58,"  was  systematic  in  that  five  parameters  characterizing 
the  body  shapes  were  perturbed  one  at  a  time  about  a  selected  "parent 
model."  Twenty-four  body  shapes  were  selected,  models  were  built,  and 
tests  were  conducted  by  towing  the  models  through  water  at  different 
speeds.  The  results  indicated  that  there  was  indeed  a  "best  shape"  for 
reducing  drag.  The  results  were  to  be  extrapolated  from  model  size  and 
speed  to  full-scale  submarine  design. 

From  the  optimization  (drag  minimization)  point  of  view,  the  one-at- 
a-time  perturbation  scheme  about  the  parent  model  represents  a  basic 
weakness  in  the  Series  58  study.  The  information  which  may  be  drawn 
from  the  study  is  the  variation  in  drag  due  to  the  variation  of  one 
parameter  while  holding  four  other  parameters  constant.  Such  information 
is  "local"  in  nature,  that  is,  the  variation  effects  are  pertinent  only 
to  the  parent  model.  Such  information  cannot  be  construed  as  "global;" 
hence,  no  claim  to  global  minimum  drag  shapes  can  be  made. 

A  design  assumption  Inferred  from  the  Series  58  study  is  that  at 
submarine  Reynolds  numbers  (Rv  -  109)  boundary  layers  are  always 
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turbulent  and  that  it  is  unrealistic  to  attempt  to  exploit  laminar 
boundary  layers  to  reduce  drag.  Without  experimental  data  it  is  diffi¬ 
cult  if  not  impossible  to  judge  the  validity  of  such  an  assumption  since 
so  many  extraneous  factors  influence  the  actual  boundary  layer  develop¬ 
ment  in  a  real  flow  situation.  The  Series  58  "best  shape"  may  be  a  low 
drag  design;  however,  it  cannot  be  inferred  from  Reference  3  that  such 
a  shape  represents  a  minimum  drag  shape  at  its  design  Reynolds  number  of 
around  109. 


1.2  Laminar  Flow  "Dolphin11  Body 

A  low-drag  shape  for  torpedo-type  Reynolds  numbers  was  reported  by 
Carmichael  in  1966  [2],  The  purpose  of  the  study  was  to  determine  if 
significant  drag  reduction  was  possible  through  shape  manipulation  alone. 
A  formal  shape-synthesizing  procedure  was  not  developed;  rather,  one 
shape  was  designed  based  on  NACA  low  drag  airfoil  data.  The  one  model, 
dubbed  the  "Dolphin,"  was  tested  extensively  by  gravity- powered  drop 
tests  in  the  Pacific  Ocean.  A  significant  drag  reduction  was  noted,  the 
"Dolphin"  having  half  the  drag  of  a  conventional  torpedo  at  similar 
Reynolds  numbers. 

The  low  drag  was  achieved  primarily  by  the  "Dolphin's"  ability  to 
maintain  a  long  run  of  laminar  boundary  layer.  This  fact  in  Itself  is 
important  since  prior  to  the  "Dolphin"  testing  it  was  generally  accepted 
that  no  significant  amount  of  laminar  boundary  layer  flow  could  be  main¬ 
tained  at  such  high  Reynolds  numbers  (Ry  =  107). 

While  the  "Dolphin"  is  not  a  minimum  drag  shape,  it  certainly 
demonstrates  the  practical  possibility  of  exploiting  laminar  boundary 
layers  and  body  shaping  in  general  to  produce  efficient  low-drag  designs, 
at  least  in  the  tested  Reynolds  number  range. 


1 .3  Outline  of  the  Text 


In  Chapter  2  the  drag  model  used  in  this  study  is  discussed  in 
detail.  The  quality  of  the  model  is  indicated  by  comparing  drag 
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predictions  with  some  of  the  experimental  data  available  in  the  liter¬ 
ature.  Chapters  3  and  4  present  the  formulation  of  the  optimization 
problem.  A  brief  discussion  of  several  possible  formulations  is  given. 
The  characteristics  of  the  drag  minimization  problem  which  lead  to  the 
selected  formulation  are  discussed.  Chapter  5  presents  the  optimum 
body  shapes  found  during  this  study.  Comparisons  with  previous  designs 
are  made,  including  comparisons  with  seme  powerful  swimmers  found  in 
nature.  Chapter  6  presents  conclusions  ana  recommendations.  The  weak¬ 
nesses  of  the  drag  model  and  optimization  procedure  are  reviewed. 
Improvements  in  both  areas  are  suggested,  along  with  proposed  future 
research. 


CHAPTER  2 


FLUID  FLOW  MODEL 


This  chapter  discusses  the  general  nature  of  the  fluid  physics  and  pre¬ 
sents  the  drag  model  as  used  in  the  optimization  problem.  Only  in 
recent  years  have  the  practical  tools  for  predicting  drag  of  arbitrary 
axisymmetric  bodies  in  incompressible  flow  appeared  in  the  literature, 
an  example  is  Reference  4.  These  methods,  while  founded  on  more  or  less 
rigorous  theoretical  considerations,  are  forced  to  rely  on  some  empiri¬ 
cal  correlations  and  simplifying  assumptions.  The  empirical  results  and 
assumptions,  along  with  the  consequential  limitations,  vary  from  one 
investigator  to  the  next;  hence,  no  one  method  of  drag  prediction  is 
regarded  as  standard.  Therefore,  it  is  the  intent  here  also  to  indicate 
the  quality  of  the  drag  model  from  which  the  quality  of  the  optimization 
results  may  be  inferred. 


General  Nature  of  the  Physical  Problem 


The  physical  phenomena  being  considered  here  are  sketched  in  Fig¬ 
ure  1.  rhe  fixed,  rigid  axisymmetric  body  is  immersed  in  an  incompressible 
f'tnd  medio  which  at  a  great  distance  from  the  body  moves  uniformly  from 
left  to  right.  The  zero  incidence  (axisymmetric)  flow  is  deflected  in 
icinity  of  the  vehicle  resulting  in  a  streamwise  pressure  gradient 
?.long  a  meridional  line  on  the  body  surface. 

The  entire  flow  field  F  is  dominated  by  inertial  (non-viscous)  forces 
except  for  a  thin  layer  next  to  the  body  surface.  Even  for  fluids  of  low 
viscosity,  such  as  water  or  air,  viscous  forces  dominate  due  to  large 
normal  velocity  gradients  which  exist  in  this  thin  boundary  layer. 
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Figure  1.  General  Incompressible  Zero  Incidence 
Flow  about  an  Axi symmetric  Body. 


Figure  2.  Simplified  Nonseparating  Flow  Considered 
in  Present  Study. 
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Because  of  its  fundamental  importance  to  vehicle  drag,  the  boundary 
layer  development  will  be  discussed  somewhat  in  detail.  Referring  to 
Figure  1,  the  boundary  layer  is  divided  into  several  streamwise  regions 
along  a  meridional  line  on  the  body  surface. 

The  laminar  region  AB  is  well  ordered  relative  to  our  own  scale 
and  can  be  readily  described  by  the  mathematics  of  continuum  mechanics. 
This  region  is  well  understood  in  the  sense  that  the  mathematical  de¬ 
scription,  with  its  assumptions,  yields  predicted  behavior  which  agrees 
well  with  that  observed  in  real  laminar  flows. 

In  the  region  BB '  the  fluid  laminas,  through  some  destabilizing 
process  not  fully  understood,  become  chaotic  on  a  macroscopic  scale. 

One  theory  [5,  6]  explains  the  process  of  transition  in  terms  of  the 
stability  of  small  disturbances  in  the  laminar  flow.  At  point  B,  down¬ 
stream  of  the  so-called  neutral  stability  point  which  is  somewhere  in 
the  laminar  region  AB,  certain  disturbance  frequencies  begin  to  amplify 
rather  than  decay.  As  the  fluid  layers  move  downstream,  the  disturbances 
amplify  and  spread  until  the  motion  of  the  entire  boundary  layer  cross- 
section  is  chaotic,  at  which  point  (B ' )  the  boundary  layer  is  said  to  be 
fully  turbulent. 

Tue  actual  transition  length  BB‘  may  be  large  or  small  depending  on 
how  fast  the  disturbances  amplify.  Some  of  the  factois  influencing  the 
location  and  extent  of  the  transition  region  for  incompressible  flow  are 
listed  below  [5]: 

1.  Free-stream  turbulence  levels  relative  to  free-stream 
mean  velocity 

2.  Streamwise  pressure  gradient  on  body 

3.  Surface  roughness  and  waviness 

4.  Body  noise  and  vibration 

The  presence  of  free-stream  turbulence,  surface  roughness  and  waviness, 
or  body  noise  and  vibration  introduces  disturbances  in  the  laminar  layer 
in  addition  to  those  already  occurring  naturally  by  self-excitation; 
the  effect  is  the  tendency  to  hasten  transition.  A  negative  pressure 
gradient  (favorable,  velocity  increasing)  tends  to  suppress  transition, 
while  a  positive  pressure  gradient  (adverse,  velocity  decreasing)  tends 
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to  hasten  traisition,  that  is,  tends  to  make  the  laminar  layer  more  un¬ 
stable.  These  qualitative  effects  have  been  verified  experimentally  [6]; 
all  quantitative  relationships  rely  on  empiricisms  to  a  greater  or  lesser 
degree. 

In  the  region  B'C  the  boundary  layer  is  fully  turbulent,  that  is, 
chaotic  on  a  macroscopic  scale.  No  rational  theory  exists  to  explain 
the  complex  fluctuating  turbulent  motion.  All  attempts  to  describe  the 
phenomena  have  relied  on  hypotheses  which  are  incomplete  without  some 
empirical  data.  One  example  is  Prandtl's  mixing-length  theory  [6];  the 
functional  form  of  the  mixing  length  is  not  established  within  the 
theory,  much  less  any  constants  within  the  functional  form. 

Turbulent  flow  is  regarded  as  the  superposition  of  mean  and  fluctu¬ 
ating  motions.  The  presence  of  the  fluctuations  has  a  decided  effect  on 
the  mean  motion,  even  when  the  fluctuating  velocities  are  a  small  per¬ 
centage  (1  -  2%)  of  the  mean  value.  Energy  is  constantly  transferred 
from  the  mean  motion  to  "large  scale"  fluctuations,  and  hence  to  pro¬ 
gressively  smaller  fluctuations;  ultimately  the  energy  is  dissipated  as 
heat.  The  apparent  effect  of  the  fluctuating  motion  is  an  increased 
resistance  by  the  mean  flow  to  fluid  deformation.  In  other  words  the 
mean  motion  behaves  as  if  the  fluid  had  an  increased  viscosity  compared 
to  the  same  fluid  in  laminar  flow.  The  apparent  viscosity  of  the  mean 
motion  depends  not  only  on  the  fluid  but  also,  apparently,  on  the  fluid 
kinematics  as  well.  The  relationship  between  the  apparent  fluid  viscosity 
and  the  fluctuating  motion  is  the  central  concept  which  defies  complete 
understanding. 

Even  with  a  lack  of  understanding  of  the  basic  mechanism  of  turbu¬ 
lent  motion,  the  mean  turbulent  flow  can  be  predicted  with  the  aid  of 
semi -empirical  relationships  [7].  The  trend  today  is  to  use  complex, 
semi -empirical  relationships,  so-called  "eddy  viscosity"  models,  to 
compute  the  apparent  increase  in  fluid  viscosity.  The  eddy  viscosity 
is  a  local  quantity  which  must  be  computed  iteratively  since  the  fluid 
kinematics  affects  the  apparent  viscosity  which  in  turn  affects  the 
kinematics.  In  all  probability,  however,  the  eddy  viscosity  concept, 
which  attempts  to  relate  apparent  viscosity  to  fluctuating  motion,  will 
be  replaced  by  mere  fundamental  turbulent  flow  shear  models  in  the  future. 
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In  the  region  CC'  the  turbulent  boundary  layer  separates  from  the 
body  surface.  By  this  it  is  meant  that  as  one  proceeds  toward  the  trail¬ 
ing  edge  along  a  meridian  on  the  body  surface,  a  point  is  reached  at 
which  the  fluid  motion  immediately  adjacent  to  the  body  reverses  direc¬ 
tion.  The  classical  view  of  laminar  separation  assumes  this  flow 
reversal  to  occur  at  a  point  in  a  steady  manner  and  to  be  synonymous 
with  zero  skin  friction.  This  view  is  not  adequate  for  turbulent 
boundary  layer  separation  which  may  be  unsteady  both  in  space  and  time 
[8].  That  is,  "spots"  of  turbulent  flow  reversal  may  occur  upstream  of 
the  fully  separated  region;  these  "spots"  may  be  unsteady  as  well.  Due 
to  a  lack  of  understanding  of  the  physics  of  turbulent  separation,  It  Is 
common  practice  to  model  turbulent  separation  as  if  it  were  a  steady 
phenomenon  occurring  at  a  point  on  the  body  meridian,  that  is,  along  an 
axisymmetric  ring  on  the  surface  of  an  axisymmetric  body. 

The  region  C'DE  is  a  fully  separated  turbulent  wake.  The  means  to 
compute  turbulent  wakes  in  the  vicinity  of  arbitrary  axisymmetric  bodies 
are  not  established,  although  wakes  behind  blunt-ended  cylinders  have 
been  predicted  successfully  [9]. 

It  should  be  apparent  that  the  complete  flow  associated  with  a  sub¬ 
merged  body  is  highly  complex.  Some  of  the  most  subtle  physics  of  n?!ure 
are  probably  associated  with  the  phenomena  of  fluid  flow.  The  incomplete 
understanding  has  forced  investigators  to  rely  to  a  greater  or  lesser 
degree  on  empirically  based  relationships,  subject  to  experimental  veri¬ 
fication,  in  order  to  predict  fluid  flow  behavior.  The  next  section 
discusses  the  simplified  flow  field,  and  its  limitations,  to  be  con¬ 
sidered  in  the  present  study  as  well  as  the  methods  used  to  predict 
such  flows. 


2.2  Simplified  Flow  Field  and  Its  Solution 

This  section  describes  the  somewhat  simplified  flow  field  and  Its 

limitations  considered  in  the  present  study.  The  methods  used  to  com¬ 
pute  the  flow  are  discussed  briefly  in  the  context  of  drag  prediction. 
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The  simplified  flow  field  considered  here  is  illustrated  in  Figure  2. 
The  flow  region  F  outside  the  boundary  layer  is  accurately  modeled  by 
invisc  *  potential  flow  theory.  The  boundary  layer  is  laminar  from  the 
forward  stagnation  point  A  to  the  transition  point  B.  Transition  is 
assumed  to  occur  at  a  single  point  B  rather  than  over  a  region.  The 
boundary  layer  is  treated  as  fully  turbulent  from  the  transition  point 
8  to  the  trailing  edge  C.  Turbulent  separation  is  assumed  not  to  occur 
so  that  the  wake  region  CDE  is  due  to  viscous  displacement  effects  only. 

In  fact,  if  a  body  has  a  separating  turbulent  boundary  layer,  as  indi¬ 
cated  by  any  suitable  separation  criterion,  then  the  relevance  of  the 
drag  model  described  here  is  not  known. 

The  flow  field  described  here  is  essentially  that  considered  by 
Cebeci  [4]  in  his  work  on  numerical  drag  prediction.  In  the  present 
study  we  have  taken  non-proprietary  versions  of  the  component  computer 
programs  developed  by  Cebeci,  A.  M.  0.  Smith,  and  their  co-workers, 
modified  them  as  required,  and  combined  them  into  a  numerical  drag 
package.  The  discussion  following  will  pertain  to  the  drag  model  as 
used  in  the  present  study  and  does  not  pertain  necessarily  to  any  drag 
package  or  component  programs  used  within  the  Douglas  Aircraft  Company. 

Inviscid  Flow  Outside  the  Boundary  Layer.  The  zero-incidence 
inviscid  flow  about  a  rather  general  axi symmetric  body  has  been  success¬ 
fully  formulated  for  numerical  solution  by  A.  M.  0.  Smith  and  his  co¬ 
workers  [10].  The  axisymmetric  problem  considered  here  is  a  subset  of 
the  more  general  geometries  and  flows  which  can  be  handled  by  Smith's 
formulation.  The  numerical  procedure  is  extremely  well  documented  [11] 
and  is  known  generally  as  the  Douglas-Neumann  method. 

In  the  context  of  the  present  problem,  Smith's  formulation  considers 
a  uniform  flow  field  in  which  is  placed  a  stationary  hypothetical  "trans¬ 
parent  image"  of  the  axisynsnetric  body.  The  body  is  temporarily  regarded 
as  "transparent"  to  the  flow  in  that  the  flow  field  remains  unperturbed. 

In  this  hypothetical  situation  the  fluid  flows  uniformly  from  left  to 
right,  entering  the  body  interior  at  some  areas  and  leaving  at  others. 
Since  the  real  body  is  impervious  to  the  fluid,  then  it  is  apparent 
that  the  hypothetical  transparent  body  must  somehow  be  made  impervious 
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also.  This  is  done  by  distributing  fluid  sources  along  the  body  surface 
so  that  at  every  point  on  the  surface  the  normal  velocity  component  is 
zero.  Stated  another  way,  the  proper  source  distribution,  which  is  . 
unique  [10],  causes  the  transparent  body  to  become  a  streamline  in  the 
flow  field.  The  superposing  of  the  uniform  flow  and  the  source  distri¬ 
bution  flow  results  in  the  inviscid  flow  about  the  impervious  body. 

The  Douglas-Neumann  method  reduces  the  problem  of  computing  a 
source  distribution  to  one  of  solving  a  linear  system.  A  modified 
Seidel  iteration  scheme  has  proven  efficient  in  solving  for  the  source 
distribution  and  hence  the  inviscid  pressure  distribution. 

The  pressure  impressed  on  the  real  body  in  a  real,  non-separating 
flow  is  the  same  as  the  inviscid  pressure  except  as  it  is  modified  by 
viscous  displacement  effects.  By  this  it  is  meant  that  the  inviscid 
streamlines  are  shifted  somewhat  both  by  the  retarded  boundary  layer  and 
viscous  wake  flows.  For  non- separating  bodies,  the  effect  of  this 
modification  to  the  inviscid  pressure  distribution  is  small.  But  it  is 
the  viscous  displacement  effect,  however  small,  that  accounts  for  the 
pressure  drag  of  a  non-separating  body. 

If  the  flow  does  not  separate  from  the  body,  then  it  is  assumed 
that  the  inviscid  pressure,  as  computed  by  the  Douglas-Neumann  method, 
is  a  reasonable  approximation  to  the  experimental  distribution.  This 
assumption  appears  to  be  justified  particularly  for  bodies  with  no 
dominating  rear  stagnation  point  in  inviscid  flow.  Examples  are  those 
with  "inflected  aftbodies,"  i.e.,  with  "semi-infinite  tailbooms"  or 
with  cusped  pointed  tails  (Figures  8  and  13).  For  these  kinds  of  bodies 
the  inviscid  flow  tends  to  free-stream  conditions  ratber  than  stagnation 

For  bodies  with  dominating  rear  stagnation  points,  which  apparently 
all  non-inflected  aftbodies  have,  the  inviscid  pressure  is  a  poor  ap¬ 
proximation  to  the  experimental  distribution  in  the  vicinity  of  the  tail, 
due  both  to  viscous  displacement  effects  and  probable  separation.  If 
the  effect  of  separation  on  the  inviscid  flow  streamlines  is  small  com¬ 
pared  to  viscous  displacement  effects,  then  it  may  be  possible  to 
approximate  the  experimental  pressure  distribution,  at  least  for  the 
purposes  of  computing  drag.  The  procedure  suggested  in  Reference  4  is 
to  replace  the  rear-stagnating  velocity  distribution  by  a  linear 
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extrapolation  of  the  inviscid  velocity  from  95%  length  to  the  trailing 
edge.  That  is,  a  straight  line  is  drawn  tangent  to  the  curve  of  in¬ 
viscid  velocity  versus  axial  distance  at  s5%  length  and  is  extended  to 
the  trailing  edge.  The  inviscid  velocity  over  the  first  95%  of  the  body 
remains  unchanged.  Using  the  modified  velocity  distribution,  the  bound¬ 
ary  layer  is  computed.  If  separation  does  not  occur,  then  the  viscous 
displacement  is  added  to  the  body  geometry.  The  inviscid  flow  about  this 
new  “viscous  body"  is  computed.  The  resulting  corrected  pressure  should 
be  a  better  approximation  to  the  experimental  distribution  than  the 
original  inviscid  flow  or  its  linearly  extrapolated  modification.  How¬ 
ever,  it  cannot  be  assumed  a  priori  that  the  corrected  pressure  is  indeed 
the  experimental  distribution  since  the  correction  depends  on  the  some¬ 
what  arbitrarily  modified  inviscid  flow.  Also,  repeated  corrections  to 
the  pressure  do  not  imply  convergence  to  the  experimental  distribution. 

A  further  comment  is  in  order  concerning  the  correction  to  the 
pressure.  The  concept  o*  viscous  displacement  is  based  on  mass  conser¬ 
vation.  That  is,  the  inviscid  flow  must  adjust  itself  to  compensate  for 
the  retarded  flow  in  the  boundary  layer  and  viscous  wake.  This  is  neces¬ 
sary  so  that  the  total  integrated  mass  flow  across  any  infinite  plane 
perpendicular  to  the  body  axis  is  conserved  when  the  flow  is  steady. 
Therefore,  the  proper  displacement  thickness  to  be  added  to  the  body 
geometry  for  the  purpose  of  correcting  the  pressure  must  be  derived 
from  mass  conservation  principles,  not  from  the  momentum  integral 
equation.  For  this  reason,  the  definition  of  displacement  thickness  6* 
given  by  equation  (2.17)  does  not  conserve  mass  flow.  The  derivation  of 
the  displacement  thickness  <5*  for  mass  conservation  in  an  external 

axisymmetric  flow  is  given  in  Appendix  A.  The  relation  between  6*  and 

aX 

6*  as  defined  by  equation  (2.17)  is 


cos  a 


+  (1  + 


2  cos  a  ,*\T 


r0 


6*): 


(2.1) 
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where  r0  and  a  are  shown  in  Figure  3.  The  quantity  5*  is  the  proper 

UA 

displacement  thickness  to  use  for  correcting  the  inviscid  pressure 
distribution.  The  value  of  6*x  is  always  less  than  that  of  6*,  partic¬ 
ularly  when  ~~  6*  is  much  larger  than  unity. 
ro 

In  summary,  the  correct  pressure  for  computing  the  boundary  layer 
and  drag  is  the  experimental  distribution.  For  inflected  aftbodies 
which  do  not  have  a  dominant  rear  stagnation  point  in  inviscid  flow,  the 
inviscid  pressure  is  assumed  to  be  a  good  approximation  to  the  experi¬ 
mental  distribution.  In  the  case  of  a  dominant  rear  stagnation  point, 
for  a  lack  of  anything  better,  the  inviscid  velocity  is  linearly  extra¬ 
polated  from  95%  length  to  the  trailing  edge  as  outlined  above.  Although 
corrections  may  improve  the  approximation  to  the  experimental  pressure 
distribution,  such  corrections  are  not  done  hero  in  order  to  conserve 
computer  time  during  optimization  runs  which  may  involve  dozens  of  drag 
evaluations. 

3oundary  Layer  Development.  Once  the  pressure  distribution  im¬ 
pressed  on  the  body  surface  is  known,  it  is  possible  then  to  proceed 
with  the  boundary  layer  computation.  A  numerical  procedure  due  to  Cebeci, 
Smith,  and  Wang  [7]  computes  planar/ axi symmetric,  laminar/ turbulent, 
incompressible/compressible  boundary  layers  using  a  variable-grid  finite- 
difference  method;  the  computer  algorithm  is  called  "Program  E7ET."  Our 
concern  here  is  the  axi symmetric,  laminar  and  turbulent,  incompressible 
problem.  The  method  of  transition  is  left  to  a  later  paragraph.  The 
method  retains  transverse  curvature  effects  which  may  be  important  for 
axi symmetric  bodies  with  thick  boundary  layers.  A  two-layer  eddy  vis¬ 
cosity  model  is  used  for  the  turbulent  boundary  layer.  Both  laminar 
and  turbulent  separation  are  based  on  the  zero  skin-friction  criterion. 
Zero  or  negative  skin  friction  is  inferred  when  the  numerical  iterative 
scheme  diverges  at  a  given  streamwise  location. 

The  coordinate  system  and  basic  notation  are  shown  in  Figure  3, 
following  Cebeci  [7].  There  is  some  redundancy  in  notation  for  two 
reasons.  In  the  past  the  notation  for  inviscid  flows  and  boundary  layer 
flows  emerged  separately.  Also,  computer  printed  output,  with  its  lack 
of  lower  case  symbols,  has  forced  some  changes  in  notation.  To  compute 


Figure  3.  Coordinate  System  and 
Basic  Notation. 


Figure  4.  Typical  Eddy  Viscosity  Distribution 
Across  Boundary  Layer. 
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inviscid  flows  and  to  describe  the  body  geometry,  the  X-Y  coordinates 
are  used;  the  reference  velocity  !)„  is  parallel  to  the  X-axis.  Boundary 
layers  are  computed  using  the  curvilinear  coordinate  system  with  x  (or  S) 
aligned  along  a  surface  meridian  and  y  normal  to  the  surface.  It  is 
convenient  also  to  use  a  radial  coordinate  r  =  r0  +  y  cos  a,  where  r0 
is  the  body  radius  and  a  is  the  angle  between  the  surface  tangent  and 
the  center  line  in  a  meridional  plane.  Ai:  any  x  location  on  the  body 
surface,  there  is  a  boundary  layer  velocity  profile  u(x,  y)  and  a 
velocity  at  the  edge  of  the  boundary  layer  ue(x).  For  incompressible 
flow  knowledge  of  the  inviscid  pressure  implies  that  the  velocity  ue(x) 
is  also  known.  Thus,  the  terms  "pressure  distribution"  and  "velocity 
distrioution"  are  used  interchangeably. 

For  incompressible  axisymmetric  flow  the  boundary  layer  equations 
[7]  may  be  written  as  follows: 

Continuity 

~  r(pu  +  p'v*’  )  +  ^  r(pv  +  pV  )  =  0  (2.2) 

Momentum 

pu|U  (pv  *  ?V  )|£  =  .&  ♦  i^r[ll|-^5TSr]  (2.3) 

where  p  is  the  fluid  density,  y  is  the  fluid  dynamic  viscosity,  p  is 
the  pressure  impressed  on  the  boundary  layer,  and  v  is  the  y-component 
of  the  boundary  layer  velocity.  Terms  containing  primes  and  overbars, 
e.g.,  p,ur,  are  time-averaged  fluctuating  quantities,  and  other  nota¬ 
tion  is  defined  above  and  in  Fiqure  3.- 

The  boundary  conditions  for  equations  (2.2)  and  (2.3)  are 

u(x,0)  ■  0  (no  slip)  (2.4a) 

v(x,0)  •  0  (no  mass  transfer)  (2.4b) 
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1 1m  u(x,y)  =  ue(x) 
y*» 


(2.4c) 


The  so-called  Reynolds  shear  stress  term  -u V  in  the  momentum 
equation  is  related  to  mean  flow  quantities  using  the  eddy  viscosity 
concept 


(2.5) 


ay 

where  e  is  the  "eddy  viscosity,"  which  is  represented  by  a  two-layer 
model  [7].  The  eddy  viscosity  for  the  inner  region,  i.e.,  the  boundary 
layer  region  near  the  wall,  is  based  on  Prandtl's  mixing-length  theory, 
so  that 

M  ■  »*  |  $  I  (2.6) 

where  £  is  the  mixing  length  defined  as 

►  « 

1  -  m  1  f-si'T  +  a*  p  >*|  <2-7> 

where  v  is  the  fluid  kinematic  viscosity,  tw  is  the  wall  shear  stress, 
and  kx  is  a  constant  equal  to  0.4.  The  outer  region  eddy  viscosity  e0 
is  given  by 

e0  =  k2  ue  6*  y  (2.8) 

which  depends  only  on  x  except  for  the  presence  of  the  intermittency 
factor  y>  where 


Y 


(2.9) 


=  [1  +  5.5(J)]'' 


and  where  6  is  the  boundary  layer  thickness  defined  as  the  y-distance 
for  which  u/ue  -  .995.  The  quantity  6*  in  equation  (2.8)  is  a  Js- 
placement  thickness  which  for  the  eddy  viscosity  model  is  defined  as 


6*  = 


(2.10) 


The  constant  k2  is  equal  to  0.0168  when  the  boundary  layer  thickness  6 
is  defined  as  the  y-distance  for  which  u/ue  =  .995.  Hence,  from  the 
wall  outward  the  eddy  viscosity  e  is  equal  to  ei  until  the  magnitude  of 
ei  equals  e0  from  which  point  outward  e  equals  e0.  A  typical  plot  of 
the  eddy  viscosity  e  across  a  boundary  layer  is  shown  in  Figure  4. 

The  purpose  for  setting  down  the  complex  expressions  for  the  eddy 
viscosity  is  to  bring  out  the  implicit  and  empirical  nature  of  the  model. 
It  is  apparent  that  the  boundary  layer  equations  (2.2)  and  (2.3)  cannot 
be  solved  until  the  eddy  viscosity  relations  (2.6),  (2.7),  (2.8),  (2.9), 
and  (2.10)  are  all  known.  It  is  also  apparent  that  equations  (2.6) 
through  (2.10)  require  the  solution  of  the  boundary  layer  before  they  can 
be  evaluated.  At  least  for  this  reason  an  iterative  scheme  is  required 
to  solve  the  turbulent  boundary  layer. 

Using  two  transformations  [7],  the  boundary  layer  equations  (2.2) 
and  (2.3)  together  with  the  boundary  conditions  (2.4),  are  non-dimen¬ 
sional  ized  and  reduced  to  an  ordinary  third-order  nonlinear  differential 
equation  with  transformed  boundary  conditions.  The  equation  is  solved 
by  an  implicit  finite-difference  method  using  variable  grid  spacing  in 
both  the  streamwise  and  normal  coordinates.  Convergence  of  solution, 
i.e.,  of  the  laminar  or  turbulent  velocity  profile,  at  each  streamwise 
station  is  based  on  convergence  of  the  transformed  wall  shear  stress. 

Transition  Prediction.  The  state-of-the-art  of  predicting  transition 
of  axi symmetric  boundary  layers  is  unclear  at  the  present  time.  An 


example  of  the  conflicting  evidence  and  a  conjectural  explanation  is 
given  in  Section  2.4.  From  the  point  of  view  of  modeling  and  simulation, 
transition  should  begin  at  the  point  where  laminar  flow  modeling  ceases 
to  be  adequate  in  some  sense.  In  the  present  study  transition  is  treated 
as  a  point  phenomenon,  that  is,  as  the  switch  which  "turns  on"  turbulence. 

The  transition  model  used  here  is  a  composite  of  a  planar-flow 
empirical  correlation  due  to  Michel  ("Michel  curve")  and  a  rather  sophis¬ 
ticated,  semi -empiri cal ,  planar-flow  correlation  ("e-to-the-nine-curve") 
due  to  Smith  and  his  co-workers  [4,  5].  The  composite  correlation,  the 
"Michel-e9  correlation,"  is  shown  in  Figure  5.  The  composite  correlation 
is  based  on  airfoil  data  taken  from  free-flight  and  low-turbulence  wind 
tunnel  tests.  The  correlation  is  between  the  momentum  thickness  Reynolds 
number  at  transition  Rq]^  and  the  running  surface-length  Reynolds  number 
at  transition  Rs]tr,  where  in  Figure  5  0  is  the  momentum  thickness,  ue 
is  the  local  velocity  at  the  edge  of  the  boundary  layer,  S  is  the  surface 
length,  and  v  is  the  fluid  kinematic  viscosity.  Also  shown  is  a  typical 
Rq  versus  Rs  plot  as  computed  by  the  boundary  layer  program  described 
in  this  study.  The  recommended  ranges  for  each  curve  are  indicated  in 
Figure  5  [4]. 

It  is  believed  that  the  axi symmetric  laminar  boundary  layer  para¬ 
meters  used  to  predict  transition  should  be  Mangier- transformed  [6]  to 
their  equivalent  planar-flow  values  before  applying  the  transition 
correlations.  However,  the  Mangier  transformation  is  explicit  only  to 
within  an  "arbitrary"  constant  which  is  not  specified  in  the  transforma¬ 
tion  itself.  The  arbitrary  constant  has  a  direct  bearing  on  the  trans¬ 
formed  boundary  layer  parameter  values  and  hence  affects  the  transition 
prediction  directly.  In  a  few  cases  that  were  checked  no  characteristic 
constant  was  found  that  consistently  improved  the  transition  prediction 
as  inferred  from  drag  data;  indeed,  in  some  cases,  the  Mangier  transfor¬ 
mation  was  grossly  detrimental  to  drag  prediction  compared  to  the  non- 
transformed  values. 

The  correlation  curve  itself  is  based  on  the  planar-flow  definition 
of  0,  which  for  incompressible  flow  is 


(2.11) 
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while  the  axi symmetric  definition  of  0  used  to  predict  transition  is 


=  f  -HJL  (1  .  JL  )  dy 

Jo  ro  ue  U  «e  1  y 


(2.12) 


where  r  is  the  radial  coordinate  of  a  point  in  the  boundary  layer  and  r0 
is  the  wall  radius.  However,  for  laminar  boundary  layers  the  value  of 
r/r0  is  nearly  1,  so  that  the  difference  in  the  definitions  (2.11)  and 
(2.12)  should  not  be  significant. 

Occasionally  the  situation  arises  in  which  laminar  separation  occurs 
before  transition  is  indicated  by  the  Michel-e9  correlation  curve.  Lami¬ 
nar  separation  is  inferred  when  the  numerical  iteration  scheme  diverges. 
In  this  situation  turbulent  reattachment  of  the  boundary  layer  is  as¬ 
sumed.  Thus,  as  far  as  the  boundary  layer  computation  is  concerned, 
laminar  separation  is  identical  to  transition. 


—  Number  and  Distribution  of  Body  Points.  It  is  appropriate  to  j 

describe  how  the  body  shape  Y(X),  which  in  general  may  exist  as  an  .] 

engineering  drawing  or  as  an  analytic  function,  is  translated  into  a  j 

table  of  points  X^,  Y(Xj),  i  =  1,  ...»  NN.  This  is  not  a  trivial  task  \ 

since  both  the  Douglas-Neumann  method  and  the  boundary  layer  program  E7ET  \ 

are  sensitive  to  the  number  and  distribution  of  body  points.  For  example, 
the  Douglas-Neumann  method  requires  a  close  spacing  of  body  points  where 
the  body  slope  or  curvature  is  changing  rapidly.  Also,  program  E7ET  \ 

requires  a  close  spacing  of  body  points  where  quantities  such  as  velocity  \ 

and  skin  friction  are  changing  rapidly.  j 

For  streamlined  bodies,  these  regions  of  rapid  change  occur  near  the  | 

leading  edge  and  in  the  vicinity  of  the  transition  point.  Since  the  \ 

transition  point  is  not  known  a  priori ,  the  points  are  closely  spaced  in  \ 

an  expanding  fashion  near  the  leading  edge  and  uniformly  spaced  else-  ; 

r 

where.  It  is  emphasized  that  the  spacing  is  established  on  the  surface  > 

meridian,  not  on  the  axis  of  symmetry.  i 

A  partial  distribution  of  points  along  a  body  surface  is  shown  in  ] 

-  —  Figure  6.  A  procedure  which  has  been  found  to  work  well  in  practice  is  1 

I 
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to  distribute  in  an  expanding  fashion  about  one-fourth  of  the  NN  points 
in  the  first  cen  to  fifteen  percent  of  the  total  surface  arc-length 
Sjotar  The  remaining  points  are  distributed  uniformly  over  the  last 

85  to  90  percent  of  Syotaj.  Near  the  leading  edge,  each  step  size  is 

equal  to  the  preceding  value  multiplied  by  the  expansion  factor  ke. 

Thus,  the  initial  step  size  is  AS0,  the  second  is  keAS0,  the  third  is 
k^ASo,  etc.  The  accumulated  arc-length  is  the  sum  of  the  three  terms, 
Ihe  initial  step  size  AS0  and  the  expansion  factor  ke  are  found  by 
solving  simultaneously  the  following  two  equations: 

AS0  (1  +  ke  +  ke2  +  ...  +  k^1-1  )  =  (FPC)  (STotal)  (2.13a) 


AS0  keN1  -  ASy 


(2.13b) 


where  N1  is  the  number  of  steps  to  be  included  in  the  first  FPC  frac¬ 
tional  part  of  SyQtal.  Equation  (2.13a)  insures  that  the  accumulated 

arc-length  of  the  first  N1  steps  is  equal  to  the  desired  fractional  part 
of  Syota-|.  Equation  (2.13b)  requires  that  the  first  uniform  step  be 

equal  to  the  previous  step  times  the  expansion  factor  ke.  The  uniform 
step  size  ASU  is  specified  by 


ASy 


sTotal 


(2.14) 


Two  sets  of  typical  values  are  given  below: 


1.  Number  of  body  points  NN  - 


FPC 

N1 

ASU 


0.10 

32/4  =  8 
•°39’3  SToUl 


32 
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AS0  =  .00387  STotal 
ke  =  1.33524 


2.  Number  of  body  points  NN  =  101 


FPC  =  0.15 

N1  =  101/4 

=  26 

ASU  =  .01134 

^Total 

AS„  =  .00299 

0 

^Total 

ke  =  1.05481 

These  two  sets  of  values  are  used  throughout  the  entire  study. 

The  set  of  NN  body  points  is  used  only  within  the  Douglas-Neumann 
method  which  computes  the  inviscid  velocities  only  at  points  midway 
between  succeeding  pairs  of  the  original  NN  body  points.  The  boundary 
layer  is  computed  only  at  the  "midpoints,"  referred  to  as  "body  sta¬ 
tions,"  by  program  E7ET.  The  body  stations  are  numbered  starting  from 
zero  so  that  the  last  body  station  number  NSTA  is  numerically  equal  to 
NN  -  2.  Throughout  this  study  the  term  "30-station"  or  "99-station" 
solution  refers  to  a  solution  obtained  using  NN  =  32  or  101  body 
points,  respectively. 

Finally,  the  algorithm  requires  60  seconds  for  30  stations 
(NN  =  32)  to  160  seconds  for  99  stations  (NN  =  101)  of  CDC  6500 
computer  time  to  compute  one  drag  value.  The  99-station  solutions  are 
more  reliable;  all  drag  values  reported  in  this  study  are  99-station 
solutions.  However,  for  an  optimization  run,  during  which  dozens  of 
drag  computations  may  be  performed,  30-station  solutions  are  used  in 
order  to  reduce  computer  costs.  The  "coarseness"  of  the  30-station 
numerical  grid  introduces  a  source  of  error  which  may  slightly  distort 
the  relative  trends  encountered  during  a  search  for  an  optimum  body 
shape.  It  is  believed,  however,  that  the  distortion  is  minor  and  that 
overall  results  are  not  affected  significantly. 

In  the  next  section  the  procedure  for  computing  drag  is  presented. 


i  t  a  iAm 


Computational  Procedure  for  Computinc 


With  the  computational  tools  outlined  in  Section  2.2,  there  are  two 
methods  immediately  available  for  computing  drag.  One  method  uses  cer¬ 
tain  flow  parameters  at  the  trailing  edge  of  the  body  in  a  formula  due 
to  Young  [12].  The  second  method  involves  integrating  the  drag  (axial) 
components  of  the  forces  acting  over  the  body  surface,  these  being 
pressure  and  skin  friction.  Both  methods  require  an  accurate  prediction 
of  the  boundary  layer  and  the  transition  location.  The  two  methods  are 
discussed  briefly  below. 


Drag  by  Young's  Formula.  Young  developed  a  formula  relating  total 
drag  to  certain  flow  parameters  at  the  trailing  edge  of  the  body.  The 
formula  in  its  non-dimensional  form  is 


where  Cq  is  the  drag  coefficient,  p  is  the  fluid  density,  U,*  is  the 

reference  velocity,  V  is  the  body  volume  so  that  V3  is  the  reference 
area,  r0  is  the  body  radius,  8  is  the  momentum  thickness  as  defined 
in  equation  (2.12),  ue  is  the  velocity  at  the  edge  of  the  boundary  layer, 
and  H  is  the  shape  factor  defined  as 

H  =  (2.16) 

where  8  is  defined  in  equation  (2.12),  and  5*  is  the  displacement  thick¬ 
ness  which  for  equation  (2.15)  is  defined  as 


S*  =  f”  L  (1  _  JL  )  (Jy 

Jo  r0  v  Ug  '  uy 


(2.17) 
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|  for  incompressible  flow.  The  definitions  of  6*.  8,  and  H  are  obtained 

s  during  the  derivation  of  the  axi symmetric  form  of  the  momentum  integral 

equation  on  which  Young's  formula  is  based.  In  equation  (2.15)  the 
subscript  T.E.  denotes  quantities  at  the  trailing  edge  of  the  body.  It 
is  noted  that  r0,  8,  ue,  and  H  are  functions  of  the  surface  length  x  so 
that  the  drag  coefficient  Cq  could  be  treated  as  Cq(x),  a  function  of  x. 
Such  a  treatment  is  not  implied  in  Young's  derivation,  although  the 
behavior  of  Cq(x)  is  of  interest.  The  derivation  o*  Young's  formula  is 
given  in  Appendix  B. 

Drag  by  Integration  of  Surface  Forces.  This  method  is  simple  in  con¬ 
cept  but  difficult  to  apply  in  practice.  Essentially,  once  the  pressure 
and  skin  friction  distributions  are  known  over  the  entire  body  surface, 
an  integration  of  the  drag  (axial)  components  should  give  the  total 
drag.  The  integral  is  given  as 

~  Drag  =  £  (p  -  pj  2irY  jjj  dX  +  Jql  tw  2ttY  dX  (2.18) 


or  in  its  non-dimensional  form  as 


where  the  first  and  second  integrals  in  (2.18)  and  (2.19)  represent  the 
contributions  due  to  pressure  and  skin  friction,  respectively.  The 
symbol  p,,,  is  tl  i-  free-stream  static  pressure,  p  is  the  local  pressure 
at  the  body  surface,  and  Cp  is  the  pressure  coefficient  defined  as 


>■¥  , 
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cp  *  j-rr  1  1  -  (  Jc  >l  (2-2°) 

b  u« 

for  incompressible  flow.  Also  Cf  is  the  skin  friction  coefficient  based 
on  the  reference  velocity  and  is  defined  as 


cf 


(2.21) 


Other  notation  is  defined  in  previous  sections. 

The  reason  equation  (2.18)  or  (2.19)  is  difficult  to  apply  In 
practice  is  that  the  pressure  over  the  body  surface  must  be  the  experi¬ 
mental  distribution,  which  Is  difficult  to  approximate  numerically  as 
discussed  in  Section  2.2.  It  appears  that  Young's  method  Is  less  sensi¬ 
tive  to  errors  in  the  pressure  distribution  than  direct  integration  of 
the  surface  forces*,  hence  Young's  formula,  equation  (2.15),  is  used  In. 
the  present  study. 

Drag  Algorithm.  The  flow  chart  for  computing  drag  In  the_  . 

present  study  is  shown  in  Figure  7.  The  input  consists  of  a  character- 

i. 

Istic  Reynolds  number  Rv  =  V5  U«/v  and  a  body  shape  Y(X).  A  geometry 
table  is  generated  using  the  procedures  associated  with  equations  (2.13) 
and  (2.14)  and  illustrated  in  Figure  6.  The  actual  numerical  values 
in  the  geometry  table  are  normalized  to  unit  length  or  unit  volume.  For 
an  optimization  run  30-station  solutions  (NN  =  32)  are  used;  for  de¬ 
tailed  drag  evaluations  99-station  solutions  (NN  *  101)  are  used. 

The  velocity  distribution  along  the  body  surface  is  computed  using 
the  Douglas-Neumann  method.  If  a  dominant  rear  stagnation  point  exists, 
then  the  inviscid  velocity  is  modified  by  linear  extrapolation  from  95 
percent  axial  length  to  the  trailing  edge. 

The  laminar  and  turbulent  boundary  layers  are  computed  by  program 
E7ET.  Transition  is  predicted  by  the  Michel-e*  correlation  or  by 
laminar  separation/turbulent  reattachment,  whichever  occurs  first. 
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Input: 


Reynolds  Number  Rv 
Body  Shape  Y(X)  v 


m  f 


Figure  7.  Flow  Chart  For  Drag  Computation. 
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If  turbulent  boundary  layer  separation  occurs,  as  indicated  by 
solution  divergence,  the  pro'’  m  aborts  and  no  drag  calculation  Is  made. 
During  an  optimization  run,  turbulent  boundary  layer  separation  Is 
treated  as  a  constraint  violation.  The  various  constraints  are  discussed 
In  Chapter  3. 

If  there  is  no  separation,  or  more  generally.  If  there  are  no  con¬ 
straint  violations  (details  In  Chapter  3),  then  the  drag  coefficient  Is 
computed  using  Young's  formula,  equation  (2.15). 

In  the  next  section  the  characteristics  of  the  drag  model  as  de¬ 
scribed  here  will  be  demonstrated  by  comparing  predictions  with  some  of 
the  data  available  In  the  literature. 


Characteristics  of  the  Drag  Model 


In  this  section  comparisons  are  made  between  predicted  and  experi¬ 
mental  drag  coefficients  using  the  drag  model  indicated  In  Figure  7. 
Transition  data  comparisons  are  also  made.  Rather  than  presenting  a 
comprehensive  study  of  the  experimental  data  available  in  the  literature. 
It  is  the  Intent  here  to  demonstrate  the  characteristics  of  the  drag 
model  used  In  the  present  study.  A  more  comprehensive  comparison  of 
predicted  and  experimental  drag  values  Is  found  in  Reference  4. 


Drag  Prediction  for  the  Laminar  Flow  "Dolphin"  Body.  A  reasonable 
approximation  to  the  "Dolphin"  Body  [2]  is  shown  In  Figure  8  along  with 
the  Inviscld  velocity  distribution.  A  typical  transition  location  Is 
indicated.  The  long,  slender  tail  boom  Is  cut  off  at  about  two-thirds 
of  the  actual  body  length.  The  experimental  data  were  obtained  from 
gravity-powered  accelerating  drop  tests  In  the  Pacific  Ocean  at  speeds 
up  to  62  knots.  Figure  9  shows  the  predicted  and  experimental  drag 
coefficients;  a  standard  torpedo  curve  [2]  Is  also  Included  for  compari¬ 
son.  The  experimental  data  has  been  corrected  to  hull  drag  values  by 
subtracting  out  the  drag  due  to  stabilizing  fins  and  the  tall  boom  [2]. 

The  agreement  between  predicted  and  experimental  values  Is  good,  with  the 
prediction  tending  to  be  optimistic.  The  role  of  the  laminar  boundary 
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layer  in  drag  reduction  is  demonstrated  both  by  prediction  and  experi¬ 
mentally  by  tripping  the  boundary  layer  at  X/L  equal  0.05.  The  drag 
coefficient  is  comparable  to  that  of  the  standard  torpedo  as  shown  in 
Figure  9.  These  data  demonstrate  that  substantial  amounts  of  laminar 
flow  are  possible  at  high  speeds  (60  knots)  in  the  ocean  environment  and 
that  the  drag  model  used  in  this  study  is  capable  of  predicting  such 
behavior  reasonably  well. 

It  is  also  of  interest  to  observe  the  trend  of  Cq  computed  at 
various  stations  along  the  body  surface.  A  typical  Cg(x)  trend  for  the 
"Dolphin"  body  is  shown  in  Figure  10.  The  drag  coefficient  Cp  increases 
monoton i cal ly  to  the  trailing  edge. 

Drag  Prediction  for  Model  4165  of  Series  58.  The  Series  58  study 
has  produced  a  recommended  best  shape  which  is  very  nearly  the  same  as 
Model  4165  of  that  series  [3].  This  body  and  its  inviscid  velocity 
distribution  are  shown  in  Figure  11.  This  body  has  a  dominant  rear 
stagnation  point  so  that  the  drag  is  computed  using  the  modified  veloc¬ 
ity  distribution  shown  in  the  same  figure;  details  of  this  modifying 
procedure  are  given  in  Section  2.3.  For  this  case  the  boundary  layer 
is  tripped  at  X/L  =  .05.  The  drag  prediction  is  for  one  Reynolds 
number,  that  corresponding  to  the  highest  test  velocity.  The  drag  co¬ 
efficient  is  computed  at  various  stations  along  the  surface  of  the  body; 
the  trend  is  shown  in  Figure  12  along  with  the  experimental  drag  value. 
As  for  the  "Dolphin"  body,  Cq  increases  monotonically  to  the  trailing 
edge. 

The  predicted  value  exceeds  the  experimental  value  by  nine  percent. 
The  predicted  value  happens  to  equal  the  experimental  value  at  the 
location  on  the  body  where  Ug/U^  is  unity.  Cebeci,  using  a  modified 
definition  of  e  in  Young's  formula,  equation  (2.15),  reports  this  be¬ 
havior  for  a  number  of  bodies  including  Model  4165  [4],  It  might  be 
inferred  that  the  drag  coefficient  is  to  be  computed  at  the  body  station 
nearest  the  trailing  edge  for  which  Ug/U,,,  is  unity.  For  the  "Dolphin," 
with  its  asymptotic  velocity,  "nearest"  would  be  interpreted  to  mean 
"at."  However,  this  idea  has  not  been  investigated;  for  the  sake  of 
consistency,  the  predicted  drag  value  will  be  that  computed  at  the  trail 
ing  edge  except  as  noted  in  the  next  paragraph. 
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Drag  Prediction  for  Murphy  Body  A2,  C4.  This  body  consists  of  a 
rounded-nose  forebody,  a  constant  diameter  midsection,  and  an  Inflected 
pointed  tail  [13].  The  boundary  layer  is  tripped  using  a  porous  strip 
about  .03L  in  width.  The  downstream  edge  of  the  strip,  located  at  X/L 
equal  0.5,  is  assumed  to  be  the  point  of  transition.  The  body  is  posi¬ 
tioned  about  0.27L  Inside  the  wind  tunnel  contraction  cone  so  that  the 
experimental  and  free-stream  inviscid  velocity  distributions  are 
different.  The  body  shape  and  its  velocity  dlstributionr  are  shown  in 
Figure  13.  Note  that  the  experimental  and  inviscid  velocity  distribu¬ 
tions  are  in  good  agreement  near  the  trailing  edge,  away  from  the 
influence  of  the  contraction  cone.  This  substantiates  to  some  degree 
the  discussion  in  Section  2.2  in  which  it  is  assumed  that  the  inviscid 
velocity  is  a  fair  approximation  to  the  experimental  velocity  distribu¬ 
tion  which  tends  to  the  free  stream  value  rather  than  rear  stagnation. 

The  drag  is  computed  at  one  Reynolds  number,  the  Cq  variation  along 
the  body  surface  is  shown  in  Figure  14  for  both  pressure  distributions 
along  with  the  experimental  range  obtained  by  wake  measurements.  The 
Cg  reaches  a  peak  value  at  about  S/S^^  equal  0.95  and  then  plunges 

rapidly.  Whe»  such  behavior  occurs,  it  is  assumed  that  the  peak  value 
is  the  proper  one  to  use,  the  rapid  plunging  apparently  indicating  a 
breakdown  in  the  method.  No  attempt  is  made  here  to  investigate  the 
underlying  reasons  for  this  behavior.  The  Douglas-Neumann  velocity 
distribution  yields  a  higher  predicted  drag  than  that  obtained  when 
using  the  experimental  distribution,  the  difference  due  mainly  to  the 
different  pressure  gradients  over  the  forebody. 

Both  predicted  values  are  in  reasonable  agreement  with  the  experi¬ 
mental  range;  the  predicted  values  are  1.8%  higher  and  5.4%  lower  than 
the  experimental  mean  value. 

Transition  Prediction  for  an  Ellipsoid.  Granville  has  reported  the 
results  of  six  different  methods  of  transition  prediction  for  axisymmetric 
bodies  using  three  experimental  pressure  distributions  on  an  ellipsoid 
with  a  fineness  ratio  (length/maximum  diameter)  of  nine  [14].  The  dif¬ 
ferent  pressure  distributions  are  obtained  by  placing  the  ellipsoid  at 
various  positions  inside  the  wind  tunnel  contraction  cone.  For  these 
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Figu'  i  13.  Murphy  Body  A2,  C4  with  Two  Velocity  Distributions. 
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tests  one  turbulence  burst  per  second  Is  taken  as  the  beginning  of  transi¬ 
tion;  the  bursts  are  detected  by  pressure  taps  in  the  body  surface. 

The  ellipsoid  and  two  of  the  three  experimental  velocity  distribu¬ 
tions  are  shown  in  Figure  15.  The  change  in  transition  location  as  a 
function  of  Reynolds  number  is  shown  in  Figure  16  for  both  velocity 
distributions.  It  is  readily  apparent  that  the  Michel-e9  correlation 
prediction  does  not  even  resemble  the  trend  of  the  experimental  data; 
in  fact,  the  correlation  predicts  nearly  the  same  transition  location 
under  all  test  conditions. 

The  apparent  contradiction  between  the  successful  "Dolphin"  drag 
predictions,  which  include  transition  prediction,  and  the  failure  of  the 
correlation  for  the  ellipsoid  has  lead  to  a  questioning  of  what  is  meant 
by  "transition"  in  the  context  of  drag  prediction. 

For  drag  prediction  the  transition  point  must  represent  the  stream- 
wise  location  on  the  body  after  which  laminar  flow  modeling  is  no  longer 
adequate.  It  is  conjectured  here  that  the  acoustical  definition  for 
transition,  i.e.,  one  turbulence  burst  per  second,  may  not  be  particular¬ 
ly  relevant  for  predicting  the  location  at  whi_n  the  boundary  layer  model 
should  "switch"  from  laminar  to  fully  turbulent.  A  recent  experimental 
study  [15]  gives  some  support  to  this  idea.  The  study  produces  a  corre¬ 
lation  among  the  average  bursting  frequency  7,  the  mainstream  velocity 
ue,  and  the  displacement  thickness  6*  using  data  from  fully  developed 
turbulent  boundary  layers  along  a  flat  wind  tunnel  wall;  the  correlation 
is  of  the  form  7  =  (constant) (ue  )/(6*).  Although  the  correlatiun 
may  not  be  directly  applicable  to  axisymmetric  boundary  layers  with 
pressure  gradients,  it  seems  reasonable  that  the  idea  of  characterizing 
a  turbulent  boundary  layer  by  using  a  correlated  bursting  frequency,  as 
opposed  to  a  fixed  frequency,  should  carry  over  to  the  axisymmetric 
case. 


Closing  Comments  about  the  Drag  Model 


The  drag  model  as  described  here  appears  to  be  reasonably  realistic. 
The  drag  predictions  may  be  either  optimistic  or  pessimistic,  apparently 


Figure  15.  Experimental  Velocity  Distributions  for  an  Ellipsoid  in  a  Contraction  Cone. 


Figure  16.  Experimental  and  Predicted 
for  Ellipsoid  Shown  in  Fig 
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depending  on  the  body  geometry  being  considered.  For  example,  the  pre¬ 
dicted  values  for  bodies  with  dominant  rear  stagnation  points  tend  to  be 
high,  but  for  inflected  aft-bodies  the  values  tend  to  be  low.  It  is 
felt  that  the  method  of  transition  prediction  tends  to  be  optimistic 
because  of  the  somewhat  optimistic  drag  predictions  of  the  "Dolphin" 
body. 

With  the  drag  model  established,  the  next  chapter  presents  the 
formulation  of  the  opiimi2ation  problem. 
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CHAPTER  3 

FORMULATION  OF  THE  OPTIMIZATION  PROBLEM 


This  chapter  contains  a  brief  discussion  of  two  possible  approaches 
to  the  optimization  problem.  The  characteristics  leading  to  the  selected 
optimization  approach  are  described;  the  constraints  are  discussed  in 
detail.  Two  optimization  methods  used  in  this  study  are  outlined. 

3.1  Functional  Optimization  —  Calculus  of  Variations 

Once  the  reference  Reynolds  number  Rv  is  fixed  for  the  zero  inci¬ 
dence  uniform  flow,  the  value  of  the  drag  coefficient  Cg  depends  on  the 
particular  shape  Y{X)  of  the  axisymmetric  body;  this  may  be  expressed  as 

CD  =  CD[Y(X)]  (3.1) 

where  X  is  the  axial  coordinate.  Equation  (3.1)  implies  that  the  drag 
coefficient  is  to  be  minimized  by  manipulating  the  function  Y(X).  Such 
a  concept  is  the  central  idea  of  the  calculus  of  variations  [18,  19]. 

The  simplicity  of  equation  (3.1)  is  somewhat  deceiving  since 
variations  must  be  performed  not  only  on  Y(X)  but  on  other  dependent 
variables  as  well.  This  is  apparent  from  equation  (2.19)  which  is 
written  here  as 


rL 


dX 


o 


(3.2) 
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where  Cp  and  Cf  are  the  pressure  and  skin  friction  coefficients  defined 
by  equations  (2.20)  and  (2.21),  respectively.  Equation  (3.2)  shows  that 

Co  is  a  function  of  Y(X),  and  the  complicated  quantities  Cp  and 

Cf  which  are  governed  by  partial  differential  equations. 

Certain  conditions  are  imposed  on  the  problem  from  the  outset. 

The  boundary  conditions  for  Y(X)  are 


Y(0)  =  0  (3.3a) 

Y(L)  =  ^Terminal  4  0  0-3b) 


where  YTenninal  is  not  necessarily  specified.  The  inequality  constraint 


Y(X)  >  0,  0  <  X  <  L 


(3.4) 


must  also  be  satisfied,  as  well  as  the  equality  constraint  tacitly  im¬ 
plied  in  equation  (3.2),  namely,  a  fixed  volume 


,L 

o 


dX 


=  V 


Specified 

Value 


(3.5) 


Furthermore,  it  may  be  necessary  to  treat  the  endpoint  X  equal  L  as  a 
variable  quantity.  A  more  thorough  analysis  will  reveal  other  con¬ 
straints  to  be  imposed  01  the  problem.  In  addition,  one  must  include 
the  physics  which  constrain  Cp  and  Cf. 

If.  the  problem  can  be  properly  formulated  using  the  integral  of 
performance  defined  in  equation  (3.2),  along  with  the  required  con¬ 
straints  and  boundary  conditions,  the  result  of  the  variational  calculus 
is  a  set  of  necessary  conditions  which  must  be  satisfied  by  the  minimum 
drag  shape.  The  set  of  necessary  conditions,  which  relate  the  dependent 
dY 

variables  Y,  Cp,  and  Cf  for  a  minimum  drag  body,  do  not  explicitly 
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define  the  optimum  body  shape  in  terms  of  the  independent  variable  X. 
However,  the  results  may  be  used  as  a  check  to  test  any  profile  which 
is  believed  to  be  a  local  optimum. 

Because  of  the  complicated  nature  of  Cp  and  Cf,  both  governed  by 
partial  differential  equations,  it  is  apparent  that  the  drag  integral, 
equation  (3.2),  cannot  be  formulated  explicitly  in  terms  of  X,  Y(X), 
and  derivatives  of  Y(X).  Without  such  a  formulation,  an  explicit  solu¬ 
tion  for  the  optimum  Y(X)  cannot  be  obtained,  at  least  not  from  a 
calculus  of  variations  analysis  alone.  For  this  reason  this  approach 
has  not  been  pursued  here. 

An  interesting  alternate  use  of  the  calculus  of  variations  has 
been  successfully  applied  to  nonseparating,  maximum  lift  airfoils  [16]. 
Rather  than  working  with  the  blade  geometry  directly,  an  optimum  pressure 
distribution  for  maximum  lift  is  obtained  from  which  the  blade  geometry 
is  uniquely  inferred.  For  axisymmetric  design,  this  approach,  called 
the  "inverse  design  problem,"  cannot  be  used  in  its  analytical  form  [17]; 
it  has  not  been  established  that  a  unique  axisymmetric  body  exists  for 
a  prescribed  pressure  distribution.  Iterative  numerical  procedures 
have  been  attempted  [17],  but  the  inverse  problem  for  axisymmetric  bodies 
does  not  appear  to  be  solved.  Hence,  this  approach  has  not  been  pursued. 


3.2  Parametric  Optimization 

If  the  calculus  of  variations  approach  were  formulated  in  terms  of 
X,  Y(X) ,  and  derivatives  of  Y(X),  and  successfully  solved,  the  solution 
for  a  fixed  Reynolds  number  Rv  would  be  an  optimum  profile  Y**(X).  For 
convenience  it  is  assumed  that  Y**(X)  is  unique.  Once  the  profile  Y**(X) 
is  known,  either  as  an  analytic  function  or  as  a  table  of  numbers,  it  is 
possible  to  approximate  it  by  a  finite  series  of  known  functions  F-j(X), 
i  =  1,  ...,  N,  each  multiplied  by  a  constant.  The  series  of  known 
functions  is  assumed  to  be  well  behaved  so  that  its  properties,  e.g., 
uniform  convergence  in  the  interval  0  <.  X  <.  L,  do  not  require  special 
consideration  here. 
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The  finite  series  approximation  may  be  expressed  as 


Y**(X) 


N 


ai**  F-j (X ) 


(3.6) 


for  OiX  <  L,  where  N  is  the  number  of  terms  in  the  series  and  the 
Fi (X) ,  i  =  1,  ....  N,  are  known  functions.  The  a^**,  1  =  1,  ...»  N, 
are  the  multiplicative  constants  which,  when  used  with  the  functions 
Fi(X),  i  *  1,  ...,  N,  yield  the  best  approximation  to  Y**(X)  in  some 
sense.  For  example,  the  a-j**,  i  =  1,  N,  in  equation  (3.6)  may 
minimize  the  error  defined  by 


(3.7) 


The  optimum  profile  Y**(X)  has  associated  with  it  a  minimum  drag  co¬ 
efficient  Cq**.  The  finite  series  on  the  right-hand  side  of  equation 
(3.6)  has  associated  with  it  a  drag  coefficient  which  can  never  be 
better  than  Cq**  since  the  finite  series  represents  a  perturbation  away 
from  the  optimum  Y**(X). 

It  is  apparent  that  by  fixing  N  and  the  F-j (X) ,  i  =  1,  ...»  N,  it 
is  possible  to  manipulate  the  multiplicative  constants  a^,  i  =  1,  ...»  N, 
also  called  "parameters,"  to  produce  an  optimum  profile 


Y*(X) 


N  * 
l  ai*  Fi(X) 
i=l 


(3.8) 


for  which  the  drag  coefficient  Cg*  is  a  minimum.  For  convenience  it 
is  assumed  here  that  the  ai*,  i  =  ...,  N,  form  a  unique  set.  It  is 
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expected  that  the  set  aj*,  i  =  1  ...»  N,  in  equation  (3.8)  is  different 

from  the  set  a^**,  i  =  1,  ....  F,  in  equation  (3.6)  since  the  former  re¬ 
sults  from  minimizing  the  drag  coefficient  Cq  while  the  latter  results 
from  minimizing  the  error  E  exemplified  by  equation  (3.7). 

Implied  in  the  above  discussion  is  the  fact  that  for  a  paramet¬ 
rically  defined  body  the  drag  coefficient  depends  on  the  number  N  of 

terms  in  the  series,  the  nature  of  the  functions  F^X),  i  =  1 . N, 

used  in  the  series,  and  the  multiplicative  constants  a-j,  i  =  1,  ...»  N. 

This  may  be  expressed  as 

CD  =  CD  [N,  a-j ,  F-j(X) ,  i  =  1,  ...,  N]  (3.9) 

This  is  true  in  particular  for  the  minimum  drag  coefficient  Cq*  asso¬ 
ciated  with  Y*(X)  in  equation  (3.8),  which  is  to  emphasize  the  fact 
that  the  minimum  drag  profile  Y*(X)  will  be  different  for  every  value 

of  N  and  every  set  of  functions  Fj(X),  i  s  1 .  N.  Therefore, 

minimum  drag  shapes  obtained  using  a  formulation  implied  by  equation 
(3.9)  can  be  regarded  only  as  the  optimum  of  a  restricted  class  of 
bodies  with  profiles  defined  by 

N 

Y(X)  =  l  a,  F,(X)  (3.10) 

i=l 

where  N  and  the  Fj(X),  i  =1,  ...,  N,  are  fixed. 

Although  the  parametric  formulation  necessarily  introduces  limita¬ 
tions  on  the  optimization  results,  it  has  been  adopted  as  the  most 
feasible  procedure  for  the  drag  minimization  problem.  Indeed,  it 
appears  that  the  functional  approach,  with  its  implicit  necessary  con¬ 
ditions,  would  require  the  use  of  an  iterative  search  strategy  of  the 
general  nature  that  we  are  considering  here  for  the  direct  problem 
solution.  A  few  comments  about  parametric  optimization  in  general  are 
given  below. 

Nearly  all  contemporary  research  in  the  area  of  optimization  methods 
is  addressed  to  the  parametric  problem  rather  than  the  functional  problem 
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of  the  classical  calculus  of  variations.  These  contemporary  methods 
are  usually  reducible  to  a  digital  computer  algorithm  so  that  they  are 
compatible  with  performance  function  models,  e.g.,  drag  models,  which 
may  already  exist  iri  digital  computer  program  form.  The  standard 
measure  of  efficiency  of  parametric  methods  is  the  number  of  perform¬ 
ance  function  evaluations  required  to  obtain  the  optimal  solution  to 
within  a  given  error  tolerance.  An  alternate  standard  is  to  compare 
the  performance  function  value  obtained  at  the  end  of  a  fixed  number  of 
evaluations. 

The  parametric  optimization  methods  may  be  broadly  classified  as 
unconstrained  or  constrained  methods.  The  generally  more  efficient 
unconstrained  methods  are  designed  to  be  used  in  a  parameter  space 
without  parameter  boundaries  (constraints).  The  generally  less  effi¬ 
cient  constrained  methods  are  designed  to  cope  with  parameter  boundaries 
(constraints)  which  divide  the  parameter  space  into  feasible  and  non- 
feasible  regions.  The  presence  of  nonfeasible  regions  may  be  due  to 
physical  considerations  or  limits  of  model  validity,  for  example.  To 
take  advantage  of  the  efficient  unconstrained  methods,  it  is  common 
practice  to  convert  a  constrained  problem  into  one  which  appears  uncon¬ 
strained  by  introduce  "penalty  functions."  A  penalty  function  arti¬ 
ficially  distorts  the  true  performance  function  "surface,"  so  that 
whenever  a  constraint  boundary  is  violated,  the  performance  function 
appears  worse  than  the  neighboring  performance  surface  in  the  feasible 
region.  The  effect  of  the  penalty  function  is  to  force  the  optimal 
solution  into  a  feasible,  hence  acceptable,  region.  It  is  generally 
desirable  to  cast  the  optimization  problem  into  one  which  is  uncon¬ 
strained  so  that  a  more  efficient  unconstrained  method  may  be  used. 

A  second  broad  classification  for  parametric  optimization  is 
gradient  versus  nongradient  methods.  This  classification  refers  to  the 
availability  of  analytical  gradients  of  the  form  3(performance  func- 
tion)/3(a-j),  where  the  aj,  i  =  1,  ...»  N,  are  the  parameters  to  be 
manipulated.  Depending  on  the  optimization  problem  this  information, 
i.e.,  the  analytical  partial  derivatives,  may  or  may  not  be  available. 
The  generally  more  efficient  gradient  methods  use  both  the  performance 
function  and  the  local  gradients  to  obtain  the  optimal  solution.  The 
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generally  less  efficient  nongradient  methods  are  designed  to  obtain  the 
optimal  solution  using  only  the  performance  function  itself.  When  the 
analytical  gradients  are  not  available,  it  is  possible  to  approximate 
the  gradients  using  finite  differences.  But  it  is  usually  more  effi¬ 
cient  to  use  a  suitable  nongradient  method  in  this  case  rather  than 
finite-difference  approximations  in  conjunction  with  an  efficient 
gradient  method  [20]. 

From  the  above  discussion  it  may  be  concluded  that  problems  which 
are  unconstrained  and  have  analytical  gradients  are  to  be  more  effi¬ 
ciently  solved  than  other  problems.  There  are  other  broad  classifica¬ 
tions  of  parametric  optimization  problems  as  well.  For  example,  the 
performance  function  and/or  the  constraints  may  be  linear  or  nonlinear 
functions  of  the  parameters.  The  constraints  may  be  expressed  as 
equalities  or  inequalities  and  may  involve  algebraic,  differential,  and 
integral  expressions.  The  constraints  may  be  explicit  or  implicit  in 
the  parameters.  The  performance  function  and/or  the  parameters  may 
be  deterministic  or  stochastic.  Tne  performance  function  and/or 
parameters  may  be  allowed  only  certain  discrete  values  rather  than 
continuously  varying  values.  The  parameter  constraint  boundaries  may 
be  convex  or  nonconvex.  The  optimal  solution  may  or  may  not  lie  on  a 
constraint  boundary.  Every  optimization  problem  will  involve  some 
combination  of  these  characteristics  and  perhaps  others  as  well. 

The  parametric  optimization  method  best  suited  to  a  particular 
problem  depends,  of  course,  on  the  characteristics  of  the  performance 
function  and  constraints  as  mentioned  above.  The  characteristics  of 
the  drag  minimization  problem  are  discussed  in  detail  in  the  next 
section.  The  discussion  leads  to  the  optimization  methods  which  are 
used  in  the  present  study. 


3.3  Characteristics  of  the  Drag  Minimization  Problem 


w 


The  implication  of  the  preceding  two  sections  is  that  the  drag 
minimization  is  to  be  cast  as  a  parametric  rather  than  a  functional 


optimization  problem.  This  section  examines  in  detail  the  characteris¬ 
tics  of  the  drag  model  and  the  constraints  in  order  to  select  appro¬ 
priate  methods  (search  strategies). 

Model  Characteristics.  From  the  parametric  optimization  point  of 
view  the  drag  model  is  a  performance  function  surface  (response  surface) 
in  an  N-dimensional  space,  where  N  is  the  number  of  independent  param¬ 
eters  (variables)  to  be  manipulated.  The  drag  model,  discussed  in 
Chapter  2,  is  essentially  a  nonlinear  numerical  "black-box"  whose  input 
is  a  set  of  parameters  a^ ,  1  =  1 ,  .. . ,  N,  and  a  Reynolds  number,  and 
whose  output  is  a  drag  coefficient  Cg.  The  parameters  imply  a  unique 
shape  Y(X)  when  N  and  F-j(X),  i  =1,  ...»  N  are  specified  in  equation 
(3.10).  Although  it  is  possible  to  treat  the  drag  prediction  in  a 
stochastic  manner  by  using  an  error  probability  distribution,  the  model 
is  treated  as  deterministic  in  the  present  study. 

The  model  is  a  "black-box"  in  the  sense  that  no  analytic  expression 
exists  relating  the  drag  coefficient  Cg  to  the  parameters  a^,  i  =1, 

....  N.  Indeed,  the  numerical  model  performs  the  same  function  as  an 
experiment,  for  example,  in  which  a  body,  whose  shape  is  Y(X)  as  implied 
by  the  aj»  i  =  1 ,  ... ,  N,  is  built  and  tested  in  a  wind  tunnel.  For 
both  the  numerical  model  and  the  hypothetical  experiment,  the  only 
information  available  for  a  given  body  at  a  fixed  Reynolds  number  is 
its  drag  coefficient  Cg.  No  analytical  gradients  3Cg/3a-j,  i  =1,  ..., 

N,  are  available  in  either  case. 

The  numerical  black-box,  with  its  lack  of  analytical  gradients,  is 
to  be  used  in  conjunction  with  a  nongradient  (direct)  search  method. 

The  procedure  of  approximating  gradients  with  finite-differences  has 
been  rejected  at  the  outset  since  it  is  believed  that  nongradient 
methods  are  more  efficient  with  "nongradient  problems"  than  finite- 
di  ■'ig  used  with  gradient  methods  [20]. 

C  •  r_  ,  int  Characteristics.  Several  statements  can  be  made  about 
the  constraints  at  the  outset.  Constraints  do  exist  for  the  drag 
minimization  problem.  The  obvious  ones  include  requirements  for  non¬ 
negative  body  dimensions,  nonseparating  flow,  and  a  fixed  Reynolds 
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number.  Less  obvious  constraints  for  a  particular  class  of  bodies 
defined  by  equation  (3.10)  include  requirements  for  a  rounded  nose  and 
no  inflection  points  on  the  forebody  (section  of  body  from  nose  to 
maximum  diameter).  These  latter  constraints  are  discussed  in  detail 
in  Chapter  4  for  two  classes  of  bodies  used  in  the  present  study. 

Other  statements  which  can  be  made  at  the  outset  are  that  the  parameters 
vary  in  a  continuous  manner  and  that  the  parameters  are  treated  as 
deterministic,  not  stochastic. 

The  procedure  for  a  performance  function  (PF)  evaluation  during 
an  optimization  run  is  to  check  for  constraint  violations  and  then  to 
compute  the  PF  if  no  violations  occur.  When  violations  do  occur,  no 
PF  evaluation  is  made;  indeed s  the  PF  value  may  not  exist  in  such  cases, 
e.g.,  negative  body  dimensions.  The  nonseparating  flow  constraint 
presents  a  special  problem;  its  violation  is  not  known  until  a  complete 
pass  is  made  through  the  drag  model.  In  an  attempt  to  reduce  computer 
time  wasted  due  to  the  occurrence  of  a  separated  boundary  layer,  two 
additional  constraints  are  checked  preceding  the  costly  boundary  layer 
computation  (75%  of  computer  time).  These  constraints  are  designed  to 
avoid  pressure  distributions  which  are  probably  conducive  to  boundary 
layer  separation.  These  constraints  are 

1.  Minimum  pressure  coefficient  Cp^  -.45 

2.  Maximum  pressure  recovery  after  the  minimum  pressure  point 
1S  CPmax  "  CPmin  £  1*0, 

Item  1  restricts  the  maximum  velocity,  which  occurs  near  the  maximum 
diameter  for  streamlined  bodies  at  zero  incidence,  to  values  less  than 
ue/Uoo  =  1.2.  Item  2  restricts  the  amount  of  deceleration  occurring 
downstream  of  the  minimum  pressure  (maximum  velocity)  point.  These 
constraints  are  not  "hard"  in  that  they  represent  reasonable  values 
but  may  rightfully  be  questioned  since  they  are  engineering  approxima¬ 
tions  to  the  separation  constraint  boundary. 

Excepting  the  Reynolds  number,  all  of  the  constraints  mentioned 
above  are  of  the  inequality  type  since  they  represent  limiting  situa¬ 
tions  or  boundaries.  The  general  conceptual  form  of  these  inequalities 
can  be  written  as 
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(3.11) 

although  the  explicit  expression  does  not  always  exist,  e.g.,  non¬ 
separating  constraint  boundary.  The  requirements  for  non-negative 
body  dimensions  are  constant  linear  constraints,  while  the  noninflected 
forebody  and  ncp.separating  flow  constraints  are  nonlinear.  Thus  the 
drag  minimization  proble...  involves  linear  and  nonlinear  inequality 
constraints;  at  '’east  one  constraint  possesses  no  explicit  form  as  given 
by  equation  (3.11). 

The  convexity  of  the  constraint  boundaries  must  be  considered. 

Regions  with  convex  boundaries,  e.g.,  interior  of  a  circle,  normally 
present  no  additional  difficulties  to  a  search  strategy.  However, 
regions  with  nonconvex  boundaries,  e.g.,  interior  of  a  cardioid,  may 
cause  a  search  strategy  to  "hang  up"  on  such  a  boundary  far  from  the 
true  feasible  optimum.  It  will  be  shown  by  graphical  means  in  Chapter  4 
that  nonconvex  boundaries  do  exist  for  at  least  one  of  the  two  classes 
of  bodies  considered  in  the  present  study. 

The  feasible  optimum  may  lie  on  the  interior  or  on  the  boundaries 
of  a  constrained  parameter  space.  Optimization  methods  which  converge 
quickly  on  the  interior  due  to  approximate  quadratic  convergence  may 
not  be  able  to  exploit  this  property  if  the  optimum  lies  on  constraint 
boundaries  (constrained  optimum).  One  reason  is  that  the  performance 
function  may  retain  dominating  first-order  properties  at  the  boundaries 
so  that  second-order  (quadratic)  characteristics  remain  insignificant. 

By  contrast,  a  performance  function  with  an  interior  optimum  will  have 
a  neighborhood  about  the  optimum  in  which  the  first-order  character¬ 
istics,  i.e.,  the  gradient  or  first  partial  derivatives,  tend  to  zero 
so  that  second-order  properties,  i.e.,  second  partial  derivatives,  tend 
to  dominate.  Since  it  is  not  known  at  the  outset  that  the  optimal 
solution  lies  on  or  off  constraint  boundaries,  it  is  appropriate  to 
consider  alternate  methods  which  work  well  in  one  situation  or  the 
ether  if  not  both. 

Reference  Reynolds  Number  —  An  Equality  Constraint.  In  previous 
discussions  the  reference  Reynolds  number  has  been  defined  as 
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Ry  =  V3  Uoo/v,  where  V  is  the  vehicle  volume  so  that  Vs  is  a  reference 
length,  Uo,  is  the  constant  vehicle  speed,  and  v  is  the  fluid  kinematic 
viscosity.  No  consideration  is  given  to  v  since  it  is  regarded  as  a 
fixed  constant  here.  A  fixed  Reynolds  number  Rv  implier  a  fixed  vehicle 
volume  V  since  'J^  is  already  specified.  Thus  Rv  is  equivalent  to  an 
equality  constraint  requiring  the  body  profile  Y(X)  to  enclose  a 
specified  volume  V.  However,  since  the  fluid  dynamics  depends  on  the 
shape  Y{X)  on d  Reynolds  number  Rv,  and  not  separately  the  volume  V, 
i.e.,  the  size  of  the  body,  and  velocity  lb  one  may  simply  scale  the 
shape  to  automatically  maintain  the  proper  volume.  In  fact,  both  V  and 
IU  may  be  scaled  so  long  as  the  proper  Rv  is  preserved.  Thus,  because 
of  the  nature  of  the  fluid  dynamics,  it  is  possible  to  exclude  consider¬ 
ation  of  the  volume  equality  constraint  from  the  optimization  strategy 
itself. 

As  mentioned  in  Chapter  1,  it  may  be  desirable  to  specify  equality 
constraints  on  quantities  other  than  volume.  For  example,  the  constant 
frontal-area  problem,  e.g.,  torpedo  design,  may  be  more  conveniently 
based  on  a  constant  maximum  diameter  Reynolds  number.  If  there  are 
neutral  buoyancy  requirements  for  the  torpedo  problem,  then  a  volume 
equality  constraint  is  still  present  and  must  be  dealt  with  by  suitable 
means,  not  necessarily  within  the  search  strategy  itself. 

Depending  on  the  application  there  may  be  other  constraints  to  be 
considered.  For  example,  submarines  must  occasionally  negotiate  chan¬ 
nels  dredged  to  a  certain  depth;  hence,  the  submarine  hull  design  is 
subject  to  a  maximum  diameter  constraint.  Another  example  is  the  design 
of  the  "lower  unit"  of  an  outboard  motor.  The  lower  unit  is  the  faired 
transmission  housing  to  which  the  propeller  is  attached  and  by  which 
the  propeller  is  powered.  The  design  of  the  lower  unit  is  subject  to 
the  constraint  of  the  space  requirements  of  the  transmission.  The 
optimization  method  should  be  able  to  cope  with  these  kinds  of  con¬ 
straints,  assuming,  of  course,  that  the  constraints  do  not  prohibit 
the  existence  of  a  feasible  solution. 

To  summarize  the  ideas  of  the  preceding  paragraphs,  the  optimiza¬ 
tion  method  used  for  the  drag  minimization  problem  must  be  capable  of 
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dealing  with  a  nonlinear  performance  function  without  the  use  of  grad¬ 
ients.  The  method  must  cope  with  nonlinear  inequality  constraints 
which  may  be  nonconvex  and  which  exist  in  implicit  or  explicit  form. 
Furthermore,  it  is  appropriate  to  consider  alternate  methods  which  are 
well  suited  to  the  optimum-on-boundary  or  optimum-on- interior  situations. 


3.4  Selected  Search  Stratf 


For  the  nongradient,  constrained,  nonlinear  drag  minimization 
problem  there  are  two  possible  approaches  using  direct  (nongradient) 
search  methods.  One  approach  is  to  use  a  method  which  operates  In  a 
nonlinear  inequality  constraint  environment.  The  second  approach  Is  to 
replace  all  inequality  constraints  with  a  suitably  constructed  penalty 
function  so  that  an  unconstrained  search  method  may  be  used.  In  either 
case  there  are  methods  which  are  regarded  as  more  efficient  than  others, 
but  the  generalization  is  not  always  valid  since  the  performance  of  a 
search  method  is  problem  dependent.  It  is  not  unusual  to  modify  the 
search  method  so  as  to  make  it  more  efficient  for  a  particular  problem. 
Some  "tailoring"  has  been  necessary  in  the  present  study;  details  are 
given  below. 

For  the  drag  minimization  problem  it  was  decided  at  the  outset  to 
select  one  promising  method  and  to  proceed  with  the  hydrodynamic  design 
problem.  Modifications  in  the  method  would  be  made  if  they  were  neces¬ 
sary  to  obtain  optimal  solutions.  No  comprehensive  experimentation 
with  various  modifications  or  different  methods  would  be  done  due  to 
the  costly  nature  of  the  performance  function  evaluation  (40  seconds  on 
the  CDC  6500).  Later  in  the  study,  however,  it  was  decided  to  try  one 
additional  method. 

Of  the  many  search  strategies  in  the  literature  for  nongradient, 
nonlinear  parametric  problems,  there  are  two  which  have  been  developed 
to  operate  in  an  environment  of  general  nonlinear  inequality  constraints 
of  the  form  given  by  equation  (3.11).  The  earlier  method  is  due  to 
Rosenbrock  (1960)  [21);  the  newer  method  is  due  to  Box  (1965)  [22]  with 
suggested  modifications  by  Guin  (1968)  [23].  In  this  study  a  slight 
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^  modification  of  the  Box-Guin  "Complex  Method"  has  been  used  extensively. 

The  method  is  described  in  a  later  section. 

The  various  nongradient  search  strategies  in  the  literature  for 
unconstrained  problems  far  outnumber  those  developed  for  the  con¬ 
strained  variety.  Examples  of  these  methods  include  an  early  approach 
due  to  Hooke  and  Jeeves  (1961)  [24],  a  method  with  approximate  quadratic 
convergence  due  to  Powell  (1964)  [25],  a  directed-hypercone  random 
search  algorithm  due  to  Wozny  and  Heydt  (1970)  [26],  and  a  recent  modi¬ 
fication  of  the  Nelder-Mead  procedure  due  to  Masters  and  Drucker  (1971) 
[27].  A  crit.'  eview  including  Powell's  method  has  been  reported  by 
Fletcher  (1965)  L28].  When  using  unconstrained  search  strategies  with 
constrained  problems,  one  common  practice  is  to  replace  the  inequality 
constraints  by  a  penalty  function  (pp.  477  -  482  of  [20]).  Of  the  many 
available  unconstrained  direct  search  methods,  Powell's  method,  be¬ 
cause  of  its  approximate  quadratic  convergence,  is  used  in  the  present 
study  in  conjunction  with  a  simple  but  general  penalty  function.  The 
details  are  left  to  a  later  section. 

Modified  Complex  Method.  A  detailed  word  flow  chart  of  the  modi¬ 
fied  Complex  Method,  as  used  in  the  present  study,  is  given  in  Appendix 
C.  The  Guin  modifications,  which  include  rtrategies  for  coping  with 
nonconvex  boundaries  and  for  generating  alternate  search  directions  when 
the  primary  direction  fails,  have  been  included  in  this  method.  A 
modification  in  the  starting  procedure  is  also  included. 

Shown  in  Figure  17  is  a  slightly  simplified  version  of  the  method; 
specifically,  the  strategy  for  nonconvex  boundaries  is  omitted  since  its 
use  has  never  been  required  during  drag  minimization  runs.  The  basic 
input  data  are  the  nu:  ->r  of  independent  parameters  N,  the  number  of 
vertices  K  in  the  complex  figure  (usually  K  =  2N),  convergence  tole¬ 
rances  e2  and  e3  (usually  e2  =  e3  =  .01),  and  the  fixed  Reynolds  number. 
For  the  initial  complex  generation,  constant  lower  and  upper  boundaries, 
aL-j  and  a^,  i  =  1,  ...»  N,  respectively,  are  also  input.  No  initial 
guess  is  needed. 

The  vertices  a_j,  j  =  1,  ...,  K,  for  the  initial  complex  figure  are 
randomly  generated  within  the  constant  rectangular  boundaries  defined 
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by  aL  ,  ay^,  i  =  1,  . ..,  N.  Each  parameter  a,  is  uniformly  distributed 
over  its  interval  ay^  -  aL^ ,  so  that  each  vertex  aj  is  uniformly  dis¬ 
tributed  over  the  enclosed  rectangular  N-dimensional  volume.  The  random 
vertex  aj  is  checked  for  feasibility;  a  nonfeasible  vertex  is  simply 
thrown  out  and  randomly  regenerated.  Whenever  a  random  vertex  is  found 
to  be  feasible,  its  performance  function  PF-j  is  evaluated.  The  process 
continues  until  a  complete  complex  is  formed. 

The  original  starting  procedure  used  by  Box  [22]  requires  an 
initial  feasible  vertex.  Each  succeeding  vertex  is  randomly  generated 
as  outlined  in  the  above  paragraph,  but  a  nonfeasible  vertex  is  moved 
halfway  toward  the  centroid  of  the  partially  completed  complex  figure. 

The  presence  of  the  initial  feasible  vortex  insures  that  the  random 
vertex  will  always  become  feasible  by  this  process.  This  procedure  is 
r  t  adequate  for  feasible  regions  which  are  highly  nonrectangular,  which 
is  the  situation  in  the  present  study.  Although  each  parameter  may  be 
allowed  large  variations,  the  actual  feasible  volume  is  a  small  fraction 
of  the  N-space  rectangular  volume  bounded  by  a_y  and  aj_.  This  effect  is 
more  prc'ounced  for  large  N.  The  effect  of  the  nonrectangularity  of  the 
feasible  region  is  to  cause  each  random  vertex  to  be  moved  half  the 
distance  to  the  centroid  many  times.  The  result  is  that  the  entire 
initial  complex  tends  to  be  clustered  in  a  relatively  small  neighborhood 
about  the  initial  feasible  vertex.  Since  the  initial  complex  is  not 
well  distributed  over  the  feasible  space,  there  is  no  global  information 
about  the  performance  function  surface.  Hence,  the  chances  of  converging 
to  the  global  feasible  optimum  are  reduced.  Furthermore.-  the  small 
scale  of  the  initial  complex  implies  small  steps  and  slow  progress  until 
the  complex  has  a  chance  to  expand.  But  the  most  detrimental  effect  is 
that  the  close  proximity  of  all  the  vertices  in  the  initial  complex 
greatly  increases  the  chances  of  premature  convergence  by  the  stopping 
criterion  used  with  this  method.  The  modified  starting  procedure 
removes  these  problems.  The  well-distributed  initial  complex  makes 
large  global  moves  at  first  and  has  a  better  chance  of  . inding  the 
global  feasible  optimum,  although  it  does  not  always  do  so,  as  results 
in  Chapter  5  demonstrate. 
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Continuing  with  Figure  17,  once  the  initial  complex  is  formed,  the 
procedure  is  to  reject  the  worst  vertex,  i.e.,  the  vertex  associated 
with  the  highest  drag  coefficient.  A  trial  vertex  is  generated  by  re¬ 
flection  through  the  centroid  of  the  nonrejected  vertices  an  amount  o 
times  the  distance  from  'he  rejected  vertex  to  the  centroid.  The 
empirical  expansion  factor  a  is  1.3  throughout  this  study,  as  recommend¬ 
ed  by  Box  [22].  If  the  trial  vertex  is  not  feasible,  it  is  repeatedly 
moved  halfway  toward  the  centroid  until  it  becomes  feasible.  Noncon¬ 
vexities  may  cause  problems  here  since  the  centroid  of  the  feasible 
vertices  may  not  lie  in  a  feasible  region.  Repeatedly  moving  a  non- 
feasible  vertex  halfway  toward  a  nonfeasible  centroid  may  prove  futile. 
Guin  suggests  at  this  point  that  the  entire  complex  be  thrown  out;  de¬ 
tails  are  in  Appendix  C. 

Once  a  feasible  trial  vertex  is  found,  its  performance  function 
PF,  i.e.,  drag  coefficient,  is  computed  and  compared  with  the  second 
worst  PF  value  of  the  complex.  The  second  worst  value  is  used  rather 
than  the  worst  to  avoid  the  situation  in  which  the  trial  PF  is  between 
the  worst  and  second  worst  values,  in  which  case  the  newly  found  vertex 
is  immediately  rejected  at  the  start  of  the  next  cycle.  This  implies 
that  the  direction  of  the  next  cycle  will  be  toward  the  point  from  which 
the  present  cycle  started.  This  situation  is  expected  when  the  complex 
is  straddling  a  local  optimum,  but  otherwise  wasted  moves  result.  If 
the  trial  PF  is  better  than  the  second  worst  PF  value,  the  trial  vertex 
replaces  the  worst  vertex;  otherwise,  the  trial  vertex  is  moved  halfway 
toward  the  centroid  and  the  new  PF  is  checked.  This  retraction  toward 
the  centroid  continues  until  the  PF  value  is  acceptable  or  until  the 
trial  vertex  enters  a  relative  e2-neighborhood  of  the  centroid.  If  the 
latter  occurs,  a  new  search  direction  is  tried  by  rejecting  the  next 
worst  vertex  and  retaining  that  previously  rejected. 

The  process  of  direction  change  continues  until  an  acceptable 
vertex  (feasible,  PF  better  than  second  worst  value)  is  found.  In  all 
cases,  the  worst  vertex  is  replaced  by  the  newly  found  vertex.  The 
stopping  condition  is  checked  and  if  it  is  not  satisfied  the  procedure 
begins  again  at  entry  point  1  in  Figure  17.  A  premature  abort  occurs 
when  a  new  acceptable  vertex  cannot  be  found.  In  such  cases,  it  is 
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—  assumed  that  an  optimum  has  been  found  even  though  the  stopping  condi¬ 

tion  is  not  satisfied.  As  with  all  search  strategies,  the  modified 
Complex  Method  converges  to  a  local  optimum;  hence,  it  may  be  necessary 
to  solve  the  problem  several  times  using  different  randomly  generated 
initial  complex  figures  in  order  to  establish  that  a  feasible  global 
optimum  has  indeed  been  found. 

Several  features  of  the  modified  Complex  Method  make  it  especially 
suited  to  the  drag  minimization  problem.  The  primary  feature  is  that 
the  method  can  cope  with  rather  general  inequality  constraints  which 
are  explicit  or  implicit  in  the  parameters.  The  method  only  requires 
a  "yes"  or  "no"  to  the  feasibility  question;  it  does  not  matter  which 
constraints  are  violated  or  how  much.  Thus  it  is  a  simple  matter  to 
cope  with  the  separation  constraint  directly.  The  complex  figure's 
ability  to  "roll"  along  boundaries  helps  to  prevent  premature  conver¬ 
gence  on  a  boundary.  Since  the  method  uses  global  features  of  the 
function  surface,  it  is  not  sensitive  to  local  irregularities  which 
might  confuse  local  gradient  methods.  The  logical  strategy  is  straight- 
forwardand  easy  to  implement  on  the  digital  computer.  The  method  as 
presented  here  is  self- starting;  if  the  optimum  is  known  to  lie  in  a 
certain  region,  that  information  can  be  exploited  at  the  outset  by 
adjusting  the  constant  boundaries  aj^  and  ay^ ,  i  =  1 ,  . . . ,  N.  No 

parameter  scaling  is  required  since  the  movements  of  the  complex  are 
automatically  scaled  to  the  range  a^  -  a^. ,  i  =  1,  ....  N. 

There  are  several  deficiencies  in  the  method  described  here.  The 
stopping  condition,  while  precisely  defined,  has  proved  economically 
costly  to  satisfy.  That  is,  many  PF  evaluations  are  required  to  es¬ 
tablish  that  the  present  best  vertex  is  a  local  optimum.  In  fact,  for 
the  drag  minimization  problem,  the  stopping  condition  has  never  been 
satisfied.  Rather,  the  search  is  aborted  after  a  large  number,  e.g., 

3N  to  4N,  of  PF  evaluations  do  not  improve  the  best  PF  value.  The  best 
vertex  is  assumed  to  be  a  reasonable  approximation  to  the  local  optimum. 
A  second  deficiency  is  the  method's  inability  to  handle  equality  con¬ 
straints;  these  must  be  handled  by  means  outside  the  search  strategy 
—  itself. 
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A  third  deficiency,  related  to  the  stopping  condition,  is  that  the 
method  lacks  quadratic  convergence  near  local  optima.  This  is  important 
if  the  local  optimum  is  interior,  not  on  a  boundary.  One  of  the  two 
classes  of  bodies  considered  in  this  study  appears  to  have  interior 
optima  (changing  with  Reynolds  number  Rv),  so  that  some  consideration 
has  been  given  to  this  weakness  of  the  Complex  Method.  It  appeal s 
appropriate  to  exploit  a  quadratically  convergent  method  for  constrained 
problems  which  have  interior  optima;  for  this  reason  Powell's  Method  of 
Conjugate  Directions  [25],  which  is  approximately  quadratically  conver¬ 
gent,  has  been  used. 
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Powell's  Method  of  Conjugate  Directions.  A  detailed  word  flow 
chart  of  Powell's  Method,  as  used  in  the  present  study,  is  given  in 
Appendix  D.  Powell's  Method  generates  search  directions  but  leaves  the 
actual  minimization  along  the  line  of  search  to  an  external  method.  A 
parabolic  interpolation  scheme  is  used  for  the  linear  minimization  in 
the  present  study;  details  are  given  in  Appendix  D. 

Shown  in  Figure  18  is  the  essential  structure  of  Powell's  Method 
of  Conjugate  Directions  as  used  in  this  study.  The  basic  input  data  are 

the  number  of  independent  parameters  N;  an  initial  feasible  guess  vector 
&0  ,  i.e. ,  a0  =  (a^,  a^,  ...»  a0fj);  1°wer  and  upper  scaling  vec¬ 
tors  a_L  and  a^,  respectively,  used  to  scale  the  search  space;  a  set  of 
linearly  independent  search  directions  £i,  £2,  ....  £N;  an  initial  step 
size  STEP  for  the  linear  search  routine  and  a  convergence  tolerance  e3. 
Mormally  the  initial  set  of  search  directions  is  the  set  of  unit  vectors 
parallel  to  the  parameter  axes.  The  scaled  parameters  vary  nominally 
between  zero  and  one. 

To  start  the  procedure  the  initial  guess  is  scaled  to  X0  ,  where 
*o  =  X°2 *  W’  using  the  re1ationshiP  xoj  =  <ao|  *  aLj)/ 

(3Ui  -  aL-j)*  The  corresponding  performance  function  value  PF0  is  eval¬ 
uated.  One  cycle  consists  of  a  linear  search  along  each  of  the  N 
direction  vectors  £2,  . ...  The  minimum  point  of  one  search  is 
the  base  point  for  the  next,  so  that  is  the  best  point  of  the  entire 
cycle.  Also,  Am  is  the  magnitude  of  the  maximum  change  in  performance 
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function  to  occur  in  any  one  linear  search.  The  corresponding  direction 
vector  of  maximum  change  is  £-jm . 

The  next  phase  of  the  method  is  to  generate  the  new  direction  vec¬ 
tor  ]i.  =  X.N  "  lo>  and  t0  compute  the  trial  point  £t=  In  +  E  Wlth  its 
performance  function  PFt.  Two  inequality  checks  determine  whether  the 
new  direction  is  promising  or  not.  If  the  new  direction  is  promising, 
an  additional  linear  search  along  it  is  made;  the  £im  direction  vector 
is  thrown  out  and  the  new  direction  vector  ji/|]i|  is  inserted  into  the 
last  position  of  the  set  of  search  directions.  Powell  has  proved  for 
quadratic  surfaces  that  this  procedure  guarantees  that  the  new  set  of 
directions  will  be  at  least  as  efficient  as  the  previous  set. 

After  deciding  to  keep  or  reject  the  new  direction,  a  convergence 
check  is  performed.  The  check  involves  the  original  base  point  iX0  and 
the  best  point  of  the  cycle.  If  convergence  is  not  achieved,  a  new 
step  size  STEP  is  computed.  The  new  STEP  magnitude  must  lie  between 
certain  reasonable  limits;  it  can  never  be  larger  than  y  (old  STEP). 

The  latter  limit  is  introduced  to  fore  the  search  to  become  more  local 
with  each  succeeding  cycle;  it  is  a  modification  of  Powell's  original 
procedure.  The  minimization  continues  cycle  by  cycle  until  the  conver¬ 
gence  criteria  are  satisfied. 

Since  Powell's  method  cannot  cope  with  constraint  boundaries,  they 
must  be  replaced  with  a  suitably  constructed  penalty  function.  The 
effect  of  the  penalty  function  must  be  to  keep  the  search  in  the  feasi¬ 
ble  region.  For  the  drag  minimization  a  penalty  function  with  the 
following  properties  is  desirable: 

1.  The  penalty  function  must  deal  with  general  explicit  or 
implicit  inequality  constraints  in  a  manner  similar  to 
the  Complex  Method. 

2.  Since  a  constraint  violation  may  render  the  drag  model 
totally  invalid,  the  penalty  function  must  generate  an 
apparent  performance  function  value  without  employing 
the  drag  model  itself.  In  other  words,  in  nonfeasible 
regions  the  response  surface  may  not  even  exist. 
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3.  The  penalty  function  should  not  distort  the  response 
surface  on  the  feasible  region. 

Item  1  implies  that  constraint  violations  are  to  be  indicated  with  a 
simple  "yes"  or  "no"  as  is  done  in  the  Complex  Method,  rather  than  indi¬ 
cating  which  constraints  are  violated  and  by  how  much.  Item  2  is  a  real 
necessity  since  some  constraint  violations,  e.g.,  non-negative  body 
dimensions,  yield  a  physically  meaningless  body  shape.  Furthermore, 
avoiding  a  drag  computation  saves  computer  time  Item  3  implies  that 
the  penalty  function  has  no  influence  on  the  feasible  region  so  that  the 
performance  function  is  the  actual  drag  value  there.  For  the  drag  mini¬ 
mization  a  performance  function  PF  which  includes  a  simple  penalty 
function  satisfying  all  three  items  above  is  defined  as 


PF 


Cq  »  no  constraints  violated 

C°NKV  +  ci|CdNKv|  »  constrai'nts  violated  (3.12) 


where  Cdnkv  is  the  nearest  known  value  of  Cp  on  the  feasible  region  and 

c1  is  a  small  positive  constant  approximately  equal  to  the  search  conver¬ 
gence  tolerance.  Equation  (3.12)  may  be  generalized  to 


PF 


PF 

actual  value 

PFNKV  +  ci |PFNKV | 
small  positive  number, 


no  constraints  violated 


constraints  violated  and 


PFNKV  /  0 


constraints  violated  and 
PFNKV  =  0 


(3.13) 


where  PF^y  is  the  nearest  known  value  of  PF  on  the  feasible  region. 
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The  constant  c1  >  0  is  arbitrary  but  is  selected  so  as  to  make  the 
ficticious  PF  surface  appear  reasonably  well  behaved. 

The  Powell  Method  used  in  conjunction  with  the  penalty  function 
implied  in  equation  (3.13)  has  several  features  useful  in  the  drag 
minimization  problem.  It  can  handle  the  general  inequality  constraints 
in  the  same  manner  a>  the  Complex  Method.  It  possesses  approximate 
quadratic  convergence  which  may  be  beneficial  for  finding  interior 
optima. 

There  are  potential  deficiencies  in  the  method.  An  initial  feas¬ 
ible  guess  is  required  so  the  method  is  not  self-starting;  for  some 
problems  locating  a  feasible  initial  guess  is  not  trivial.  The  param¬ 
eter  space  should  be  scaled  so  that  each  parameter  has  an  "equal  inter¬ 
est"  in  the  performance  function.  The  method  uses  local  information  in 
its  moves  so  that  local  irregularities  may  confuse  the  search  strategy. 
However,  one  motivation  for  using  a  "local"  method  in  the  drag  minimiza¬ 
tion  is  to  study  the  migration  of  local  minima  with  Reynolds  number. 

The  stopping  condition  as  indicated  in  Figure  18  is  as  costly  to  satisfy 
as  that  for  the  Complex  Method.  From  a  practical  point  of  view  it  may 
be  adequate  to  test  for  performance  function  convergence  but  not  param¬ 
eter  convergence. 


3.5  Properties  of  the  Optimal  Solution 

After  obtaining  an  optimal  solution,  four  properties  of  that  solu¬ 
tion  must  be  considered.  These  include  uniqueness,  a  global  ''ersus 
nonglobal  solution,  and  the  sensitivity  of  the  optimum  to  off-design 
conditions.  The  fourth  property  emerges  when  the  optimal  solution  is 
obtained  using  a  finite  search  and  a  finite  stopping  condition;  it  is 
the  proximity  of  the  reported  solution  to  the  true  local  minimum.  These 
properties  are  of  concern  in  all  optimization  problems,  but  they  will  be 
commented  on  below  in  terms  of  a  numerical  "black-box"  performance 
function  to  which  a  finite  search  has  been  applied.  This  is,  of  course, 
the  situation  for  the  drag  minimization  problem. 
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The  t e/'i.  :  .que"  is  being  used  in  an  engineering  sense.  Optimal 
solutions  w.  :ch  are  spatially  far  apart  in  an  N-parameter  space  but 
wnich  are  close  in  their  performance  function  values,  say  to  within 
the  convergence  tolerance,  are  regarded  as  nonunique  solutions.  The 
tern  "global"  refers  to  the  best  of  all  the  feasible  local  optima,  but 
it  has  meaning  only  after  uniqueness  is  established.  For  "black-box" 
models  there  are  no  rigorous  procedures  for  establishing  whether  an 
optimal  solution  is  either  unique  or  global.  This  is  true  for  optima 
on  the  interior  or  on  a  '  ;undary.  At  best  these  properties  can  only  be 
indicated  by  solving  the  same  optimization  problem  several  times  using 
different  starting  conditions.  This  procedure  has  been  used  in  the 
present  study. 

The  question  of  sensitivity  has  real  practical  importance.  Es¬ 
sentially,  it  is  desirable  to  know  how  quickly  the  optimal  solution 
degrades  in  performance  at  off-design  conditions.  Such  conditions  occur 
due  to  variations  in  the  parameters,  variations  in  assumed  fixed  con¬ 
ditions,  e.g.,  Reynolds  number,  and  variations  in  the  model,  e.g., 
transition  prediction.  These  off-design  conditions  are  examined  in 
this  study,  although  not  uniformly  for  every  optimum  body  design. 

The  proximity  of  the  reported  solution  to  the  true  local  optimum 
can  be  interpreted  in  two  ways.  From  the  design  point  of  view  proxim¬ 
ity  of  performance  is  emphasized;  from  the  optimization  point  of  view 
spatial  proximity  as  well  as  performance  proximity  are  important.  The 
latter  statement  is  true  because  it  is  of  interest  to  know  how  efficient 
a  search  strategy  is  in  seeking  out  local  minima.  When  using  finite 
search  strategies  with  "black-box"  models,  the  local  optimum  is  usually 
never  known  exactly.  One  procedure  is  to  fit  an  analytic  quadratic 
surface  locally  in  the  neighborhood  of  the  best  solution  [20].  The 
ninimi  1  of  this  analytic  surface  is  found  accurately;  the  true  minimum 
is  known  to  within  a  finte  but  -.taller  tolerance.  A  side  bene  it  of 
the  quadratic  surface  fit  is  the  immediate  approximate  information  *f 
local  curvature  behavior  useful  in  sensitivity  studies. 

To  work  reliably  the  above  procedure  requires  a  reasonable  distri¬ 
bution  of  experiments  (point? )  in  a  ne  ghborhood  about  the  suspected 
local  optimum.  In  an  N-parameter  space  at  least  N  +  N(N  +  l)/2 
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experiments  are  required.  Usually  the  requirement  of  a  "reasonable 
distribution"  forces  the  generation  of  some  additional  points  not  ob¬ 
tained  during  the  actual  optimization  run  so  that  an  extra  cost  is 
incurred  using  this  procedure. 

A  less  thorough  and  less  costly  approach  for  determining  the 
proximity  of  the  best  point  to  the  true  local  optimum  to  examine 
several  points  randomly  distributed  on  the  surface  of  an  N-hypersphere 
whose  center  is  the  reported  best  point.  The  radius  of  the  normalized 
hypersphere  may  be  the  normalized  convergence  tolerance  for  example. 

The  procedure  Is  to  randomly  generate  a  direction,  orthogonal  to  any 
previous  directions;  and  to  test  the  performance  function  value  at 
both  ends  of  the  hypers pheri cal  diameter  parallel  to  the  generated 
direction.  The  random  direction  is  rejected  and  regenerated  if  both 
ends  are  in  nonfeasible  regions.  If  no  better  point  is  found  after 
two  or  three  random  directions  have  been  checked,  then  the  confidence 
in  the  assumption  that  the  reported  optimum  is  near  the  true  local 
optimum  has  increased.  However,  if  a  better  point  is  found,  then  the 
radius  is  immediately  increased  and  the  procedure  is  repeated.  An 
alternate  procedure  is  to  re-center  the  same  hypersphere  on  the  new 
best  point  and  to  repeat  the  procedure.  If  a  bettt;  point  is  found 
on  the  second  hypersphere,  then  it  is  assumed  that  the  reported  optimal 
solution  is  not  a  particularly  good  approximation  to  the  true  local 
optimum.  Eut  the  point  may  still  be  acceptable  if  the  performance 
function  value  differences  arc.  small.  In  the  present  study  proximity 
checks  are  made  using  this  procedure  rather  than  the  quadratic  surface 
fit. 


3.6  Comment  on  Optimization  Philosophy 


The  optimization  p-  ’rr.ophy  emerging  in  this  chapter  is  summarized 
in  this  section.  Essenl.ally  the  drag  minimization  is  to  be  solved 
through  the  interactions  of  two  digital  computer  programs,  each  a 
"black-box"  to  the  other.  The  drag  model  and  search  strategies  are 
independent  of  each  other;  hence,  independent  improvements  can  be  made 


in  one  without  affecting  the  other.  This  idea  has  been  established 
from  the  outset  so  that  as  better  drag  models  or  search  strategies 
become  available,  they  may  be  incorporated  into  the  drag  minimization 
package  with  only  minor  programming  changes.  So  as  to  not  overempha¬ 
size  the  "black-box"  relationship,  experience  has  shown  that  it  is 
usually  beneficial  to  tailor  the  search  strategy  somewhat  to  the  partic¬ 
ular  problem  in  order  to  make  the  complete  package  more  efficient. 

The  Complex  Method  and  Powell's  Method  with  the  penalty  function 
given  by  equation  (3.13)  represent  diverse  search  strategies.  The 
diversity  should  lend  a  degree  of  confidence  to  the  determination  of 
uniqueness  and  global  optimality.  The  Complex  Method,  with  its  ran¬ 
domly  and  globally  distributed  initial  complex  figure,  has  some  chance 
of  finding  the  global  optimum,  if  there  is  only  one.  Powell's  Method, 
with  its  local  movements,  should  be  able  to  "track"  a  local  minimum 
drag  shape  with  changing  Reynolds  number.  Two  other  benefits  of 
Powell's  Method  for  interior  optima  include  efficient  convergence  and 
the  approximate  knowledge  of  the  local  curvature  at  the  optimum  point. 

There  are  two  points  of  view  regarding  the  results  of  the  drag 
minimization  studies.  One  is  the  design  point  of  view  in  which  the 
emphasis  is  on  the  performance  of  the  design  and  its  sensitivity  to 
off-design  conditions.  The  other  is  the  optimization  point  of  view 
which  includes  consideration  of  the  above  and  also  uniqueness,  globality, 
proximity  to  the  true  optimum,  and  search  efficiency.  Both  points  of 
view  are  retained  in  the  results  to  follow  but  not  uniformly  with  every 
minimum  drag  shape. 
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CHAPTER  4 

PARAMETRIC  BODY  PROFILES 


This  chapter  defines  the  two  classes  of  bodies  used  in  the  present 
study.  Certain  inequality  constraint  expressions  are  also  presented. 

A  parametric  definition  of  body  shapes  more  general  than  that  given  by 
equation  (3.10)  may  he  expressed  as 


Y(X) 


Y 


l  aiFi(X),  G(X) 


(4.1) 


where  N  is  the  number  of  independent  parameters  a-j  associated  with  the 
known  function  F^(X).  The  known  function  G(X)  is  present  to  satisfy 
certain  boundary  conditions  built  into  Y(X).  Equation  (4.1)  implies 
that  Y(X)  may  not  be  a  simple  linear  combination  of  known  functions. 

The  procedures  used  to  derive  the  expressions  for  Y(X),  i.e., 
equation  (4.1),  are  essentially  those  reported  by  Granville  (1969)  [29]. 
The  idea  is  to  divide  the  body  into  sections  each  of  which  is  described 
by  a  low  degree  polynomial.  From  hydrodynamic  v  nsideratio  s  the  com¬ 
plete  body  is  to  be  continuous  through  second  derivatives;  for  example, 
a  discontinuity  in  curvature  can  cause  a  "pressure  spike"  (local  region 
of  highly  accelerated  flow)  to  occur.  The  low  degree  polynomial  of 
c'ch  section  is  completely  specified  in  terms  of  its  boundary  condi¬ 
tions,  some  of  which  are  fixed  and  some  of  which  are  free  to  be  manipu¬ 
lated.  Those  which  may  be  manipulated  are  the  parameters  to  be  varied 
during  an  optimization  run.  Furthermore,  for  each  body  section  only 
two  boundary  conditions  are  free  so  that  two-dimensional  constraint 
boundaries  can  be  plotted.  The  complete  derivations  are  left  to  Ap¬ 
pendix  E;  th  results  are  reported  here. 
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4.1  Five- Parameter  Rounded-Nose,  Pointed-Tail  Body  [29] 

An  example  of  this  body  is  shown  in  Figurel9.  The  forebody  (OiX  < 
Xm)  is  described  by  a  fourth-degree  polynomial,  the  aftbody  (Xm  <  X  i  L) 
by  a  fifth-degree  polynomial. 


Figure  19.  Rounded-Nose  Pointed  Tail  Body. 


The  six  dimensional  parameters  shown  in  the  figure  are  listed  below: 

1.  Rn  =  Radius  of  curvature  at  nose 

=  l/(d2X/dY2)  at  X  =  0 

2.  D  =  Maximum  diameter 

-  Y(Xm) 

3.  Xm  =  Axial  location  of  maximum  diameter  D 
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4.  Kj  =  Curvature  at  Xm 

-  d2Y(Xm)/dX2 

5.  St  =  Profile  slope  at  tail 

=  dY(L)/dX 

6.  L  =  Overall  body  length 

The  rounded-nose,  i.e.,  infinite  slope,  is  a  built-in  boundary 
condition,  although  a  zero  radius  of  curvature,  Rn  =  0,  is  allowed.  The 
pointed-tail,  i.e.,  finite  slope,  is  also  a  built-in  boundary  condition. 

The  six  parameters  listed  above  can  be  reduced  to  five  which  are 
nondimensional ,  hence  the  "five-parameter"  designation.  These  are 
listed  below. 

1.  rn  =  Nondimensional  radius  of  curvature  at  nose 


=  C4Xm/D2]  Rn  =  [4xmf£]  Rn/L  (4.2a) 

2.  fr  =  Fineness  ratio 

=  L/D  (4.2h) 

3.  xm  =  Nondimensional  axial  location  of  maximum  diameter  D 

=  Xm/L  (4.2c) 

4.  kj  =  Nondimensional  curvature  at  Xm 

=  [-2XJ/D]  Kx  =  [-2xmfr]  L  (4.2d) 

5.  Sj-  =  Nondimensional  slope  at  tail 

=  [-2(L-Xm)/0]  St  =  C-2(l-xm)fr]  St  (4.2e) 


All  nondimensional  parameters  are  defined  so  that  they  are  normally 
non-negative.  The  particular  nondimensional izing  expressions  emerge 


; 


during  the  derivations  given  in  Appendix  E,  The  length  L  is  not  a  free 
parameter  to  be  manipulated  since  it  must  be  scaled  to  satisfy  the 
fixed  Reynolds  number. 

The  analytical  expressions  for  this  five-parameter  body  are  given 
below: 


1 .  For  0  i  X  <  Xra  (forebody): 


i 


11*1 . 

2^:  CrnF1  (x)  +  kjFjjfc)  +  G(x)]7 

(4.3) 

where  x 

•  X/Xm 

(4.4a) 

F,(x) 

=  -2x(x  -  l)3 

(4.4b) 

F.(x) 

=  -x2 (x  -  l)2 

(4.4c) 

G(x) 

=  x2(3x2  -  8x  +  6) 

(4.4d) 

2.  For  Xm  <  X  i  L  (pointed  aftbody): 


^  ^  [4  F,(*>  +  (^f1)2  k,  F2(x)  +  S(x)]2  (4.5) 

where  x  =  (L-X)/(L-Xm)  (4.6a) 

Fx(x)  =  -x2(x  -  l)3  (4.6b) 

F2(x)  =  -x3(x  -  l)2  (4.6c) 

G(x)  =  x3(6x2  -  15x  +  10)  (4.6d) 


The  constraint  boundaries  imposed  on  the  five  nondimensional  param¬ 
eters  are  listed  below.  Some  are  obviously  "pre-judgments"  influenced 
by  previous  hydrodynamic  experience. 


■A^*L****  C_ 
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The  complicated  constraints  are  listed  last  and  are  considered  in 
detail . 

1.  Non-negative  radius  of  curvature  Rn  at  nose,  or 

rn^0  (4.7a) 

2.  Nonpositive  curvature  Kt  at  maximum  diameter,  or 

^>0  (4.7b) 

3.  Real  location  for  maximum  diameter  so  that 

0  <  xm  <  1  (4.7c) 

4.  Reasonable  fineness  ratios,  i.e.,  not  conducive  to 
separation,  so  that 

fr  >  some  positive  constant,  say  2.5  (4.7d) 

5.  Nonpositive  slope  St  at  tail,  or 

st  >  0  (4.7e) 

6.  No  inflection  points  on  forebody 

7.  No  or  only  one  inflection  point  on  aftbody 


The  low  degree  polynomials  for  the  two  body  sections  together  with 
constraints  1  through  7  above  imply  that  the  body  profile  is  always 
non-negative. 

Noninflected  Forebody.  The  forebody,  which  is  described  by  a 
fourth-degree  polynomial,  may  have  zero,  one,  or  two  inflection  points. 
The  analysis  [29]  given  in  Appendix  E  leads  to  a  set  of  two  simultaneous 
nonlinear  algebraic  equations  which  must  be  solved  by  numerical  itera¬ 
tion.  The  solution  gives  the  rn  versus  kj  curve  along  which  one  limiting 
inflection  occurs  on  the  forebody  in  that  the  second  derivative  touches 
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zero  but  does  not  change  sign.  The  convex  rn  versus  boundary  is  shown 
in  Figure  20;  the  non-negative  boundaries  are  due  to  equations  (4.7a) 
and  (4.7b)  above.  The  existence  of  two  inflections  in  the  labelled 
region  is  demonstrated  graphically  [29]. 


Figure  20.  Feasible  Region  for  Noninflected 
Rounded -Nose  Forebody  Section. 


Noninflected  or  Inflected  Pointed  Aftbody.  The  aft.body,  which  is 
described  by  a  fifth-degree  polynomial,  may  have  zero,  one,  two,  or 
three  inflection  points.  The  analysis  [29]  given  in  Appendix  E  leads 
to  a  set  of  two  simultaneous  nonlinear  algebraic  equations  which  must 
be  solved  by  numerical  iteration.  A  singularity  leads  to  another  set 
of  equations  which  can  be  solved  directly.  The  resulting  s|  versus 
[(1  -  xm)/xiiJ2ki  curves  are  shown  in  Figure  21.  A  nonconvexity  arises 
when  both  noninflected  and  singly  inflected  aftbodies  are  allowed.  The 
curved  boundary  represents  the  limiting  inflection  in  that  the  second 
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derivative  touches  zero  but  does  not  change  sign.  The  number  of  inflec¬ 
tions  in  each  region  has  been  established  by  Granville  [29]. 


Figure  21.  Feasible  Region  for  Noninfiected  or  Singly 
Inflected  Pointed  Aftbody  Section. 


4.2  Eight-Parameter  Rounded-Nose,  Tailboom  Bodj 


An  example  of  this  body  is  shown  in  Figure  22.  The  forebody  (0  < 
X  <.  Xm)  is  described  by  a  fourth-degree  polynomial  as  in  Section  4.1 
above,  the  midbody  (Xm  <  X  i  X-j)  bj  a  fifth-degree  polynomial,  and  the 
tailboom  aftbody  (Xj  <  X  <.  L)  by  a  fifth-degree  polynomial. 


*rt*s  •> 
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Figure  22.  Rounded-Nose  Tailboom  Body. 


The  nine 

dimensional  parameters  shown  in  the  figure  are  listed  below: 

1. 

Rn 

= 

Radius  of  curvature  at  nose 

= 

l/[d2X/dY2]  at  X  =  0 

2. 

D 

= 

Maximum  diameter 

= 

Y(Xm) 

3. 

xm 

r 

Axial  location  of  maximum  diameter  D 

4. 

= 

Curvature  at  Xm 

= 

d2Y(Xm)/dX2 

5. 

*i 

= 

Axial  location  of  inflection  point 

77 

6.  R*  =  Profile  radius  at  Xj 

=  Y(Xi) 

7.  S.j  =  Profile  slope  at  X< 

=  dY(Xi)/dX 

8.  T  =  Terminal  profile  radius 

»  m 

9.  L  =  Overall  body  length 

The  rounded-nose,  i.e.,  infinite  slope,  is  a  built-in  boundary  condition 
as  in  Section  4.1  above,  and  Rn  =  0  is  allowed.  The  tailboom  has  built- 
in  boundary  conditions  of  zero  slope  and  curvature,  and  T  =  0  is  mathe¬ 
matically  allowable.  The  zero  slope  and  curvature  at  X  equal  L  implies 
that  a  cylindrical  tailboom  extension  may  be  added  withov*  loss  of 
profile  continuity  through  the  second  derivative. 

The  nine  parameters  listed  above  can  be  reduced  to  eight  which  are 
nondimensional ,  hence  the  "eight-parameter"  designation.  These  are 
listed  below: 

1.  rn  =  Nondimensional  radius  of  curvature  at  nose 

=  (4X[n/D2)  Rn  =  ( ^xmf r )  Rn/L  (4.8a) 

2.  fr  =  Fineness  ratio 

=  L/D  (4-8b) 

3.  xm  =  Nondimensional  axial  location  of  maximum  diameter  D 

=  Xm/L  (4.8c) 

4.  kj  =  Nondimensional  curvature  at  Xm 

=  (-2X>)  Kj  =  (-2x*fr)  K,  L 


(4.8d) 
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5.  x-j  =  Nondimensicnal  axial  location  of  inflection  point 


=  Xi/L  (4.8e) 

6.  Tj  =  Nondimensional  profile  radius  at  X-j 

=  (2/D)  Ri  =  (2fr)  Ri/L  (4.8f ) 

7.  s-j  =  Nondimensional  profile  slope  at  X-j 

=  C-(Xi  -  Xm)/(D/2  -  Ri)]Si 

=  [-2fr(xj  -  xm ) / ( 1  "  ^ i )HS-j  (4.8g) 

8.  t  =  Nondimensional  terminal  profile  radius 

-  (2/D)  T  =  (2fr)  T/L  (4.8h) 


All  nondimensional  parameters  are  defined  so  that  they  are  normally  non¬ 
negative.  The  particular  nondimensional izing  expressions  emerge  during 
the  derivations  given  in  Appendix  E.  The  length  L  is  not  a  free  param¬ 
eter  to  be  manipulated  since  it  must  be  scaled  to  satisfy  the  fixed 
Reynolds  number. 

The  analytical  expressions  for  this  eight-parameter  body  are  given 
below.  The  expressions  for  the  forebody  (0  <  X  <  Xm)  are  identical  to 
those  given  by  equations  (4.3)  and  (4.4)  in  Section  4.1  above;  they  are 
not  repeated  here. 

1.  For  Xm  <.  X  <.  X-j  (midbody): 
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where 

x  = 

(X j  -  X)/^  -  xm) 

(4.10a) 

Ft(x)  = 

-|x3(x  -  l)2 

(4.10b) 

F2(x)  = 

x  -  x3(3x2  -  8x  +  6) 

(4.10c) 

G(x)  = 

x3(6x2  -  15x  +  10) 

(4.10d) 

For  X- 

i  <  X  £  L 

(tailboom  aftbody): 

Y(X) 

L 

=  fl 
2fr  l1 

5  ,  t  n  r  /..x  ,  (l-nMl-Xi] 

+  (n  -  »  r-(xi  +  I^JTnT 

-  Si  F2(x) 

4 

(4.11) 

where 

x  = 

(L  -  X)/(L  -  X,) 

(4.12a) 

Fj  (x)  = 

1  -  x3(6x2  -  15.  +  10) 

(4.12b) 

F2(x)  = 

-x3(3x2  -  7x  +  4) 

(4.12c) 

The  constraint  boundaries  imposed  on  the  eight  i  ndimensional  param¬ 
eters  are  listed  below.  Because  the  forebody  is  iden- ical  to  that  for 
the  five- parameter  body,  some  of  the  constraints  listed  here  ara  dupli¬ 
cates  of  those  found  in  Section  4.1. 

1  Non-negative  radius  of  r-  *ature  Rn  at  nose,  or 

rn  2.  0  (4.13a) 

2.  Nonpositive  curvature  Kx  at  maximum  diameter,  or 

kx>0  (4.13b) 

3.  Compatible  locations  of  maximum  diameter  and  inflection 
point  so  that 


(4.13c) 
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4.  Reasonable  fineness  ratios,  i.e.,  not  conducive  to 
separation,  so  that 

fr  >.  some  positive  constant,  say  2.5  (4.13d) 

5.  Compatible  radii  at  maximum  diameter,  inflection 
point,  and  trailing  edge,  so  that 

0  <  t  <  r-j  <  1  (4.13e) 

6.  Nonpositive  profile  slope  S-j  at  inflection  point,  or 

s,-  >  0  (4.1 3f ) 

7.  No  inflection  points  on  body  except  at  and  L 

The  low  degree  polynomials  for  the  three  body  sections  together  with 
constraints  1  through  7  above  imply  that  the  body  profile  is  always 
non-negative.  The  inflection  point  requirement,  item  7,  implies  that 
all  three  body  sections  are  noninflected  except  at  Xj  and  L.  The 
analysis  for  the  noninflected  forebody  is  Identical  to  that  for  the 
five-parameter  body  so  that  Figure  21  applies. 

Noninflected  Midbody.  The  midbody,  which  is  described  by  a  fifth- 
degree  polynomial,  may  have  zero,  one,  two,  or  three  inflection  points, 
but  no  more  than  two  on  the  interval  Xm  <  X  <  X-j  since  one  is  fixed  at 
X  equal  X^.  The  analysis  in  Appendix  E  leads  to  a  boundary  curve  de¬ 
fined  by  two  explicit  parametric  equations.  A  singularity  leads  to  an 
additional  equation.  The  resulting  convex  Si  versus  [ ( x -j /xm  -  l)2/ 

(1  -  r-j )3kj  curves  are  shown  in  Figure  23.  The  number  of  inflections 
in  each  region  is  established  by  a  general  analytical  demonstration. 

Noninflected  Tail  boom  Aftbody.  The  tail  boom  aftbody,  which  is 
described  by  a  fifth-degree  polynomial ,  may  have  zero,  one,  two,  or 
three  inflection  points,  but  no  more  than  one  on  the  interval  X-j  <  X  < 

L  since  one  is  fixed  at  X  equal  Xj  and  another  at  X  equal  L.  The 
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Figure  23.  Feasible  Region  for  Noninflected  Midbody  Section 
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Figure  24.  Feasible  Region  for  Noninflected 
Tailboom  Aftbody  Section. 
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analysis  in  Appendix  E  gives  the  linear  convex  [(1  -  r-j ) (1  -  x-j)/ 

(x-j  -  xm)(ri)]sj  versus  t/r-j  boundaries  shown  in  Figure  24.  The  number 
of  inflections  in  each  region  is  established  by  a  general  analytical 
treatment. 


4.3  Closing  Comment 

It  is  apparent  that  the  two  classes  of  bodies  considered  in  this 
study  are  constrained  in  advance  to  be  well  behaved  according  to  previ¬ 
ous  hydrodynamic  experience.  The  profiles  are  continuous  through  all 
derivatives  except  at  a  finite  number  of  points  which  join  body  sec¬ 
tions;  at  such  points  the  profiles  are  continuous  through  second 
derivatives.  The  discontinuous  third  derivative  at  these  points  implies 
that  the  curvature,  while  continuous,  may  change  rapidly.  Such  behavior 
causes  local  regions  of  accelerated  flow.  It  is  interesting  to  note 
that  the  search  strategies  described  in  Chapter  3  have  exploited  this 
phenomenon  to  help  minimize  drag,  as  shown  in  the  next  chapter. 
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CHAPTER  5 

RESULTS  AND  COMPARISONS 


This  chapter  presents  results  obtained  using  the  drag  model  and 
search  strategies  described  in  the  preceding  chapters.  Some  attention 
is  also  given  to  the  hydrodynamic  performance  of  powerful  swimmers 
found  in  nature.  Overall  conclusions  are  reserved  for  Chapter  6. 


5.1  Body  D-54  and  the  "Dolphin" 

The  impressive  performance  of  the  "Dolphin"  [2],  discussed  in 
Chapters  1  and  2,  represents  the  standard  of  -omparison  for  the  eight- 
parameter  tail  boom  body.  It  is  of  interest  to  know  whether  or  not  a 
body  superior  (lower  Cq)  to  the  "Dolphin"  can  be  found.  Using  the 
Complex  Method,  an  optimization  run  made  early  in  the  study  has  pro¬ 
duced  a  body  with  a  drag  coefficient  Cq  about  25  percent  lower  than 
that  of  the  "Dolphin"  at  similar  Reynolds  numbers.  The  resulting  body, 
called  "D-54,"  and  its  velocity  distribution  are  shown  in  Figure  25 
along  with  the  "Dolphin"  profile.  Body  D-54  is  the  54th  function  (Cq) 
evaluation  of  the  optimization  run. 

Body  D-54  is  characterized  by  a  long  run  of  laminar  boundary  layer 
flow  over  the  forward  two-thirds  of  the  body.  The  small  velocity  gra¬ 
dients  over  the  forebody  help  to  reduce  skin  friction;  however,  the 
same  near-zero  gradients  have  a  neutrally  stabilizing  effect  on  the 
laminar  boundary  layer,  i.e.,  the  absence  of  an  accelerating  boundary 
layer  increases  the  chances  of  early  transition.  At  the  midsectinn  the 
laminar  boundary  layer  is  approaching  conditions  for  transition  as  in¬ 
dicated  by  the  Rq  versus  Rs  trajectory  in  Figure  26.  Transition  is 
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suppressed,  however,  by  !ocally  accelerating  the  flow;  this  is  accom¬ 
plished  by  allowing  the  body  curvature  to  change  rapidly.  After 
transition,  which  is  indicated  by  laminar  separation/turbulent  re¬ 
attachment,  the  turbulent  boundary  layer  survives  the  run  of  adverse 
velocity  gradient  and  enters  the  terminal  accelerating  region.  Such  a 
region  helps  to  suppress  boundary  layer  separation  as  suggested  by  the 
skin  friction  distribution  plotted  in  Figure  27;  we  are  equating  nonzero 
skin  friction  and  nonseparating  flow. 

Perturbation  studies  on  body  D-54  have  been  made.  The  procedure 
is  to  randomly  generate  orthonormal  direction  vectors  6  in  Euclidean 
8-space.  For  example,  all  components  of  61  are  randomly  generated; 
each  succeeding  vector  has  one  less  random  component  so  that  the  re¬ 
maining  components  may  be  used  to  satisfy  orthogonality  requirements. 
Three  orthonormal  vectors  used  in  this  study  are  shown  in  Table  1.  The 
reported  optimum  body  shape  is  represented  by  the  set  of  parameters  £*, 
where  for  the  eight-parameter  tailboom  body 

a*  =  (a*,  a*,  a*,  a*,  a*,  a*,  a*,  aj) 

=  (f*.  xj,  k*.  r*,  r*,  st,  x*.  t*)  (5.1) 

Perturbations  a/  about  the  reported  optimum  a*  are  generated  by 

aj  =  aj  (1  ±  R6j),  j  -  1 . 8  (5.2) 

where  R  is  the  magnitude  of  the  perturbation  and  6j,  j  =  1,  ....  8,  are 
the  perturbation  direction  components.  Body  D-54  and  six  3%  perturba¬ 
tions  are  given  in  Table  2  using  the  directions  of  Table  1  and  R  =  .03. 

The  perturbed  drag  coefficient  values,  which  include  two  6X  per¬ 
turbations  not  included  in  Table  2,  are  shown  in  Figure  28;  it  is 
obvious  that  body  D-54  is  suboptimal ,  a  result  of  prematurely  stopping 
the  search  after  65  function  evaluations.  Since  the  normalized  gradients 
3CD/3(perturbation  direction)  are  on  the  order  of  one,  significant  im¬ 
provements  in  the  minimum  Cq  should  be  possible. 


Friction  Coefficient 


Table  1.  Three  Random  Orthonormal  Directions  in  Euclidean  8-Space. 


Table  2.  Body  D-54  with  Six  3%  Perturbations. 
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Hqure  28.  Perturbation  Results  for  Body  D-54. 
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5.2  All-Turbulent  Body  1-36  and  the  Series  58  Stud\ 


The  philosophy  of  the  Series  58  study  [3],  discussed  in  Chapters  1 
and  2,  is  to  design  low  drag  bodies  on  the  assumption  that  laminar  flow 
cannot  exist  in  the  operating  environment  at  submarine- type  Reynolds 
numbers.  No  attempt  is  made  here  to  judge  the  validity  of  this  assump¬ 
tion;  rather,  it  is  of  interest  to  know  whether  or  not  an  all -turbulent 
body  can  be  found  with  a  Cq  lower  than  that  of  the  best  of  the  Series 
58  bodies. 

Using  the  Complex  Method  and  the  five-parameter  pointed  tail  body, 
three  optimization  runs  with  different  initial  complex  figures  have 
been  made.  The  boundary  layer  is  tripped  at  X/L  equals  .05;  for  all 
three  runs  the  Reynolds  number  is  fixed  at  Ry  equals  5  x  10s  which  is 
the  upper  end  of  the  Series  58  test  Ry  values.  Two  of  the  three  runs 
converged  after  27  function  evaluations;  the  third  run  was  near  con¬ 
vergence  but  was  aborted  prematurely  after  38  function  evaluations  due 
to  reasons  external  to  the  search  strategy. 

The  significant  result  of  the  three  runs  is  the  fact  that  the 
response  surface  for  all-turbulent  bodies  is  quite  flat.  That  is,  for 
fairly  wide  variations  in  the  parameters,  the  drag  coefficient  varies 
little.  This  is  evident  in  Figure  29  which  shows  nine  parameter  sets 
plotted  in  the  parameter  space.  All  parameter  sets  have  corresponding 
Cq  values  within  one  percent  of  the  best  design.  It  is  evident  that  low 
drag  all-turbulent  bodies  are  not  critically  dependent  on  shape.  Hence, 
if  laminar  boundary  layers  cannot  be  exploited  to  minimize  drag,  then 
means  other  than  profile  shaping  must  be  used  to  reduce  drag,  e.g., 
polymer  injection. 

The  best  design,  body  "1-36,"  and  its  velocity  distribution  are 
shown  in  Figure  30  along  with  the  profile  of  model  4165,  the  best  of 
the  Series  58  study.  The  Cq  value  for  body  1-36  is  0.020,  which  is 
about  the  same  as  that  of  model  4165  according  to  the  drag  model  used 
in  this  study  The  streamwise  velocity  gradient  is  small  over  most  of 
the  body  length  to  reduce  skin  friction. 
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5.3  A  Series  of  Laminar  Bodies  at  Three  Reynolds  Numbers 

Of  primary  interest  is  the  change  of  minimum  drag  shape  over  a 
wide  range  of  Reynolds  numbers.  An  underlying  topic  of  fundamental 
importance  is  the  uniqueness  of  a  minimum  drag  shape  at  its  design 
Reynolds  number.  In  this  section  we  report  the  results  of  optimization 
runs  made  at  three  Reynolds  numbers  using  the  five- parameter  pointed- 
tail  body  described  in  Chapter  4.  The  selected  Reynolds  numbers  are 
Ry  =  5  x  10s ,  1.6  x  107,  and  5  x  107  which  correspond  to  nominal  vehicle 
volumes  of  1.1,  37,  and  1130  cubic  feet  traveling  at  a  speed  of  35  knots 
in  water. 


Low  Ry  Body  6-35.  Using  the  Complex  Method  a  low  drag  body  shape 

has  been  obtained  at  Ry  =  5  x  106.  The  optimization  run  terminates 
prematurely  after  43  function  evaluations  due  to  reasons  external  to  the 
search  strategy.  The  best  body  shape  occurs  on  the  35th  function  eval¬ 
uation  of  the  run;  its  predicted  Cg  value  is  0.0054  at  the  design  Ry. 

Body  6-35  and  its  velocity  distribution  are  shown  in  Figure  31. 

The  body  is  quite  streamlined  and  has  a  long  run  of  laminar  boundary 
layer  flow  over  the  forward  three-quarters  of  th^  body.  The  forebody  is 
essentially  of  the  Reichardt  type  since  the  veldfcity  gradient  is  nearly 
zero  over  most  of  the  forward  half  of  the  body.  Iri  a  manner  similar  to 
the  tail  boom  body  D-54,  the  flow  is  locally  accelerated  starting  at 
about  X/L  equal  0.6  to  suppress  transition,  which  is  predicted  by  lami¬ 
nar  separation/turbulent  reattachment.  The  effect  is  readily  seen  in 
the  Rq  versus  R$  trajectory  plotted  in  Figure  32.  The  skin  friction 
distribution  plotted  in  Figure  33  suggests  that  the  turbulent  boundary 
layer  is  on  the  verge  of  separating  near  the  trailing  edge. 

Perturbation  studies  have  been  made  to  determine  if  body  6-35  is 
near  a  local  minimum.  The  parameters  of  the  reported  optimum  in  a 
Euclidean  5-space  are 


a*  =  (a*,  a*,  a*,  a*,  a*) 

■  <•  *  t-  H2> 


(5.3) 


"  ^  V".- 


I  R0  =  Log  6  ue/v 


7 


Log  Rs  =  Log  S  ug/v 


Figure  32.  RQ  versus  R$  for  Body  G-35 
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As  for  the  eight-parameter  body  0-54,  random  orthonormal  directions, 
shown  in  Table  3,  are  used  to  generate  perturbations 

aj  =  aj^  *  R6j>*  j  =  •••»  5  (5.4) 

where  R  is  the  perturbation  magnitude  and  the  6j,  j  =  1,  ...»  5,  are 
the  components  of  the  perturbation  direction.  The  parameters  for  body 
6-35  and  six  Z%  perturbations  (R  =  .03)  using  the  directions  in  Table  3 
are  shown  in  Table  4.  The  perturbation  results,  shown  in  Figure  34, 
reveal  Immediately  that  body  G-35  is  suboptimal.  The  normalized  gradi¬ 
ents  3(Cg)/3(perturbation  direction)  are  on  the  order  of  one  so  that 
significant  improvement  in  the  minimum  Cg  should  be  possible. 

Midrange  Ry  Body  H-62.  Using  the  Complex  Method  a  low  drag  body 

shape  has  been  obtained  at  Ry  =  1.6  x  107.  The  optimization  run  termi¬ 
nates  after  80  function  evaluations,  the  last  improvement  occurring 
on  the  62nd  evaluation. 

Body  H-62  and  its  velocity  distribution  are  shown  In  Figure  35;  the 
Cg  value  is  .0059  at  the  design  Ry  value.  The  body  is  somewhat  "fatter" 
and  more  pointed  than  the  low  Ry  body  G-35.  A  long  run  of  laminar 
boundary  layer  flow  is  maintained  by  continuously  but  mildly  accelerat¬ 
ing  the  flow  over  the  forward  two-thirds  of  the  body.  Transition  is 
predicted  by  the  Michel-e9  correlation,  as  indicated  by  the  R0  versus  R$ 
trajectory  plotted  in  Figure  36.  The  skin  friction  distribution  shown 
in  Figure  37  suggests  that  the  turbulent  boundary  layer  is  on  the  verge 
of  separating  near  the  trailing  edge  in  a  manner  similar  to  body  G-35. 

Perturbation  studies  have  been  made  on  body  H-62  using  the  pertur¬ 
bation  directions  given  in  Table  3.  The  procedure  is  the  same  as  for 
body  G-35.  The  parameters  for  body  H-62  and  six  3%  perturbations 
(R  =  .03)  are  shown  in  Table  5.  The  perturbation  results  are  shown  in 
Figure  38;  direction  6^  reveals  chat  body  H-62  is  suboptimal  and  that 
improvement  in  the  minimum  Cg  value  is  possible. 

High  Ry  Bodies  F-57  and  F2-49.  Using  both  the  Complex  Method  and 
Powell's  Method  two  distinct  low  drag  shapes  have  been  obtained  at 
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Table  3.  Three  Random  Orthonormal 
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Directions  in  Euclidean  5-Space. 
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Ry  =  5  x  107.  One  optimization  run  using  the  Complex  Method  is  termi¬ 
nated  after  90  function  evaluations,  the  best  design  occurring  on  the 
57th  iteration  (body  F-57 ) .  A  second  run  using  the  Complex  Method  is 
terminated  after  70  function  evaluations  producing  a  different  best 
shape  on  the  49th  iteration  (body  F2-49).  A  third  run  using  Powell's 
Method  is  terminated  after  75  function  evaluations  with  a  best  shape 
quite  similar  to  body  F2-49.  None  of  the  three  runs  converges  according 
to  the  formal  definitions  associated  with  the  methods. 

Body  F-57  and  its  velocity  distribution  are  given  in  Figure  39; 
the  drag  coefficient  is  .0076.  This  "fat"  body  has  some  hydrodynamic 
similarity  to  body  H-62  in  that  it  is  pointed  and  has  a  long  run  of 
accelerated  laminar  flow  over  the  forward  two-fifths  of  the  body.  Un¬ 
like  the  lower  Ry  body  shapes,  there  is  no  dominant  locally  accelerated 
flow  to  suppress  transition.  The  absence  of  the  effect  is  seen  in 
Figure  40  in  which  the  Re  versus  R$  trajectory  approaches  the  Michel -e9 
correlation  curve  in  a  monotonic  manner,  ultimately  crossing  the  curve 
to  predict  transition.  The  skin  friction  distribution  shown  in  Figure 
41  suggests  that  the  turbulent  boundary  layer  is  on  the  verge  of  sepa¬ 
rating  near  the  trailing  edge. 

Body  F2-49  and  its  velocity  distribution  are  given  in  Figure  42; 
the  drag  coefficient  is  .0073  so  that  it  is  superior  to  body  F-57  on 
the  basis  of  minimum  Cq.  This  shape  also  has  hydrodynamic  similarity 
to  body  H-62  in  that  a  small  run  of  nearly  constant  velocity  flow 
precedes  a  large  region  of  accelerated  flow  which  helps  to  suppress 
transition.  This  effect  is  seen  in  Figure  43  in  which  the  Re  versus 
R$  trajectory  approaches,  veers  away,  and  approaches  again  the  Michel- 
e9  correlation  curve;  its  ultimate  crossing  predicts  transition.  In  a 
manner  similar  to  all  the  preceding  pointed  tail  bodies,  the  skin 
friction  distribution  shown  in  Figure  44  suggests  that  the  turbulent 
boundary  layer  is  on  the  verge  of  separating  near  the  trailing  edge. 

It  is  interesting  co  note  that  each  of  the  two  low  drag  shapes 
exploits  a  different  feature  of  the  midrange  Ry  body  H-62  to  minimize 
drag.  That  is,  the  pointed-nose  body  F-57  maintains  laminar  flow  by 
continuously  accelerating  the  flow  over  the  forebody;  the  rounded-nose 
body  F2-49  suppresses  transition  by  locally  accelerating  the  flow 


ty  Ratio 


y  F2-49  with  Velocity  Distribution. 


,1  =  X/L 


Laminar 
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starting  at  about  X/L  equal  0.2.  Even  though  the  hydrodynamic  means 
are  somewhat  different,  transition  occurs  at  about  the  same  axial  loca¬ 
tion  for  both  bodies. 

The  existence  of  these  two  distinct  low  drag  shapes  brings  the 
question  of  uniqueness  into  the  foreground.  In  an  attempt  to  establish 
uniqueness,  i.e.,  the  existence  of  a  finite  number  of  distinct  low  drag 
shapes,  a  third  optimization  run  has  been  made  using  Powell's  Method; 
the  idea  is  to  see  whether  a  different  search  strategy  converges  to  an 
existing  solution  or  produces  yet  another  low  drag  design.  If  the 
former  occurs,  then  confidence  in  the  uniqueness  of  the  solutions  has 
increased.  If  the  latter  occurs,  then  nonuniqueness  appears  likely  and 
further  testing  is  required. 

The  results  of  the  Powell  run  are  best  seen  by  observing  the  over¬ 
all  parameter  migrations  of  all  three  optimization  runs  as  shown  in 
Figure  45.  The  data  shown  are  the  initial  and  finai  designs  for  each 
run;  the  arrowhead  indicates  direction  of  overall  movement.  For  the 
runs  using  the  Complex  Method,  the  initial  design  shown  is  the  best 
vertex  of  the  initial  complex  figure.  It  can  be  seen  that  both  fir^l 
designs  are  similar  in  their  parameter  values  except  for  the  radius  of 
curvature  at  the  nose  rn.  It  should  be  noted  that  small  differences  in 
xm  are  exaggerated  on  the  ka[(l  -  xm)/xm]2  axis.  The  result  of  the 
Powell  ser.rch  is  quite  close  t.o  body  F2-49  and  appears  to  be  converging 
on  the  rn  -  kx  boundary.  The  conclusion  to  be  drawn  is  that  the  low 
drag  designs  are  distinct  local  minima,  finite  in  number,  and  hence 
unique. 

Variation  of  Cq  with  Reynolds  Number.  One  sensitivity  analysis  of 

interest  to  the  hydrodynamic  designer  is  the  variation  of  drag  coeffi¬ 
cient  with  Reynolds  number.  The  variation  of  Cq  over  a  wide  range  of 
Ry  is  shown  in  Figure  46  for  low,  midrange,  and  high  Ry  bodies.  At  its 
own  design  point  each  low  drag  body  has  the  lowest  Cq  of  the  three 
values.  The  designs  are  not  sensitive  to  reasonably  wide  variations  In 
Ry,  at  least  according  to  predictions  of  the  drag  model  used  in  this 
study. 
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5.4  Hydrodynamic  Performance  of  Powerful  Swimmers 


Data  on  the  swimming  speeds  of  dolphins  (porpoises)  have  appeared 
in  the  literature  since  the  mid  1930* s .  Until  the  early  1960's  all 
data  indicated  that  dolphins  either  1)  produced  several  times  the  power 
believed  possible  by  their  musculature,  or  2)  were  able  to  reduce  their 
flow  resistance  by  unknown  means  to  values  several  times  lower  than 
those  of  similar  man-made  devices  [31,  32],  However,  recent  studies 
[33,  34]  have  shown  through  reasonably  well  controlled  experiments  that 
the  dolphin  does  not  have  low  resistance  to  flow;  indeed,  typical  drag 
coefficients  are  about  the  same  as  that  of  a  torpedo  with  an  all-turbu¬ 
lent  boundary  layer.  Stated  another  way,  there  are  rigid  laminar  flow 
bodies  with  drag  coefficients  about  half  that  of  the  dolphin  animal. 

A  tested  example  of  such  a  man-made  device  is  the  "Dolphin"  [2]  dis¬ 
cussed  in  Chapters  1  and  2. 

Table  6  compares  drag  data  taken  from  the  references  indicated  for 
a  nominal  Ry  of  5  x  10s.  Data  for  the  porpoise  species  are  usually 
reported  on  a  dimensional  basis  so  that  direct  comparisons  are  difficult 
to  make,  item  1  of  Table  6  being  a  rare  exception.  It  is  apparent  from 
the  table  that  the  Stenella  attenuata  has  a  drag  coefficient  more  like 
torpedoes  and  all -turbulent  bodies  than  like  laminar  flow  devices.  It 
must  be  concluded  that  porpoises,  as  hydrodynamic  performers,  are. 
mediocre. 

There  are  at  least  three  reasons  why  the  earlier  conclusions  con¬ 
cerning  the  hydrodynamic  abilities  of  the  dolphin  are  erroneous  [33,  34]. 
Unusually  high  speeds  of  dolphins  swimming  near  moving  ships  may  be 
explained  by  assisted  locomotion  since  it  is  known  that  dolphins  are  able 
to  derive  thrust  from  the  moving  fluid  near  the  ship.  Some  speeds  have 
been  deduced  from  dolphin  sightings  by  observers  on  moving  ships  over  a 
quarter-mile  distance  [32].  However,  the  combination  of  ship  motion  and 
wave  motion  creates  an  illusion  of  fast  speed  which  has  been  estimated 
to  cause  errors  of  25%  [34].  A  more  subtle  source  of  error  is  the 
duration  of  time  over  which  the  speeds  are  recorded.  Dolphins  are 
capable  of  high  power  output  and  hence  high  speeds  for  short  periods  of 
time,  perhaps  5  to  10  seconds,  during  which  the  muscle  tissue  goes  into 


Table  6. 


Porpoise  Drag  Data  Compared  with  Rigid 
Devices  for  Nominal  Ry  of  Five  Million 


Animal 
or  Device 
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1.  Porpoise,  Stenella 
attenuata 


.00239 

to 

.00401 


2.  Typical  Torpedo 


.0196 


Drag  due  to  appendages 
subtracted  out.  Data 
from  coasting  tests 
(1966)  [34]. 


[2] 


3.  Model  4165 
of  Series  58 


.00255 


All  turbulent  boundary 
layer.  Tow  tank  data  [3], 


4.  Laminar  Flow 
"Dolphin" 


.00109  .0092  Data  from  gravity-pow- 

to  to  ere.,  drop  tests  [2J. 

.00164  .0138 


*  A  =  body  wetted  surface  area 
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oxygen  debt.  One  series  of  tests  [33]  showed  that  the  dolphin  attained 
a  speed  of  16.1  knots  for  7.5  seconds  but  only  6  knots  on  a  continuous 
basis.  For  this  example  the  ratio  of  peak  power  to  continuous  power 
is  approximately  (16.1/6)3  =  19.2.  Apparently  the  concept  of  peak 
power  has  not  been  considered,  or  at  least  has  not  been  stated,  in 
earlier  work. 

In  comparing  the  shapes  of  powerful  swimmers  to  those  of  rigid  low 
drag  bodies,  one  important  fact  must  be  emphasized:  the  optimum  shape 
for  swimming  -.limals  will  be  different  from  the  optimum  shape  of  low 
drag  man-made  devices.  The  underlying  reason  is  that  the  animals  and 
devices  are  optimally  shaped  with  respect  to  different  performance 
functions.  The  animal  is  shaped  so  as  to  maximize  ics  chances  for 
survival,  which  is  a  rather  complex  performance  function  indeed!  On 
the  other  hand,  the  shapes  of  similar  man-made  devices,  as  proposed  to 
date,  are  optimal  with  respect  to  the  relatively  uncomplicated  perform¬ 
ance  function  of  drag.  The  powerful  swimmer  is  shaped  so  as  to  attempt 
to  simultaneously  minimize  drag,  maximize  propulsion,  maximize  energy 
reserves,  minimize  the  effect  of  interfering  objects  such  as  eyes, 
mouth,  and  gills,  and  minimize  its  metabolism  rate  while  swimming 
quickly.  To  extremize  any  one  of  these  items  without  regard  for  the 
others  would  certainly  change  the  shape  of  the  animal. 

Of  equal  importance  is  the  fact  that  the  powerful  swimmers  are 
optimally  shaped  subject  to  a  different  set  of  constraints  than  those 
for  man-made  devices.  Animals  are  constrained  to  utilize  propulsion 
mechanisms  which  oscillate  rather  than  rotate  since  all  parts  of  the 
animal  must  be  connected  with  blood  vessels  and  nerves.  There  is  also 
••  the  subtle  constraint  of  maintaining  a  favorable  surface-area-to-volume 
ratio  so  that  the  organism  may  function  properly  [35],  Man-made  de¬ 
vices  are  practically  constrained  to  use  rotating  mechanisms  for  pro¬ 
pulsion  since  oscillating  mechanisms  of  rigid  parts  apparently  cannot 
be  made  to  operate  efficiently.  Implied  here  is  the  additional  practi¬ 
cal  constraint  of  rigid  construction  and  structural  integrity. 

The  point  to  be  made  here  is  that  it  is  not  particularly  meaningful 
to  compare  the  shapes  of  powerful  swimmers  and  low  drag  rigid  devices. 
Each  is  the  result  of  attempting  to  extremize  different  performance 
functions  subject  to  different  constraints. 
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CHAPTER  6 

CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions  and  recommendations  given  in  this  chapter  are  based 
on  results  reported  in  Chapter  5  and  Appendix  F.  At  the  outset  it  must 
be  said  that  significant  vehicle  drag  reduction  is  possible  through  shape 
manipulation.  The  present  method  has  produced  low  drag  bodies  with  drag 
coefficients  one-fourth  to  one-third  below  that  of  the  low  drag  "Dolphin" 
at  similar  Reynolds  numbers.  Of  equal  significance  is  the  fact  that  the 
minimum  drag  shape  is  a  strong  function  of  Reynolds  number.  For  the 
five-parameter  body  over  a  one  order-of-magnitude  range  of  Reynolds 
number  (5  x  106  <  Ry  <  5  x  107)  the  corresponding  optimum  fineness  ratio 
(L/D)  ranges  nominally  from  8.5  to  3.5.  The  corresponding  location  of 
maximum  diameter  ranges  nominally  from  Xm/L  =  .75  to  .45. 

Over  the  reported  Ry  range  all  optimum  body  shapes  exploit  laminar 
boundary  layers  to  reduce  drag.  The  experimental  evidence  of  the  "Dol¬ 
phin"  [2]  demonstrates  that  substantial  laminar  flow  does  exist  in  the 
ocean  environment  at  speeds  up  to  60  knots  (Ry  above  107).  It  must  be 
concluded  that  laminar  flow  is  a  practical  means  for  reducing  drag  at 
these  Reynolds  numbers  and  that  proper  body  shaping  can  use  laminar 
flow  effectively. 

If  laminar  flow  is  prevented  due  to  extraneous  factors,  such  as 
body  roughness,  then  it  appears  that  the  body  shape  is  not  particularly 
critical  in  reducing  drag  so  long  as  it  is  reasonably  streamlined.  This 
conclusion  is  based  on  the  all-turbulent  drag  minimization  study  at 
Ry  =  5  x  10s  in  which  nothing  better  than  the  best  of  the  Series  58 
study  [3]  is  obtained.  Stated  another  way,  if  laminar  flow  cannot  be 
exploited  to  reduce  drag,  then  further  drag  reduction  below  present 
design  values  must  be  accomplished  by  means  other  than  profile  shaping 
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alone.  One  alternate  means  is  polymer  injection  in  the  boundary  layer. 

Numerous  optimization  runs  made  at  Ry  =  5  x  TO7  using  two  search 
strategies  produced  a  number  of  low  drag  shapes.  Fine  tuning  by  hand 
reveals  that  there  is  apparently  a  unique  global  minimum  drag  shape. 

The  global  minimum  has  high  sensitivity  to  early  transition;  hence,  sub- 
optimal  solutions  without  this  sensitivity  are  more  desirable  from  the 
hydrodynamic  design  point  of  view.  Alternatively,  additional  constraints 
may  be  imposed  on  the  problem  to  avoid  such  undesirable  characteristics. 

The  modified  Complex  Method  used  in  this  study  has  performed  well. 

It  operates  in  a  constrained  environment  without  difficulty.  Since  it 
moves  on  global  information,  the  method  can  cope  with  errors  in  the 
performance  function  which  do  not  obliterate  global  trends.  The  method 
is  well  suited  for  design  problems  without  critical  convergence  tole¬ 
rances.  For  the  drag  minimization  problem,  one  may  expect  to  use  ION  to 
15N  minutes  of  CbC  6500  computer  time  to  obtain  a  minimum  drag  body 
using  the  recommended  30-station  solutions,  where  N  is  the  number  of 
parameters. 

Powell's  Method  along  with  the  penalty  function  used  in  this  study 
can  only  be  used  effectively  with  99-station  (3  minute)  solutions  since 
the  normal  30-station  (1  minute)  solutions  introduce  enough  error  to 
confuse  this  locally  moving  method.  The  method’s  primary  use  for  the 
present  problem  is  fine  tuning;  one  may  expect  to  use  9N  to  12N  minutes 
of  CDC  5500  computer  time  for  each  cycle  of  fine  tuning  using  the  recom¬ 
mended  99-station  solutions,  where  N  is  the  number  of  parameters. 
Normally,  fine  tuning  by  this  method  is  not  necessary;  when  it  is  re¬ 
quired,  one  cycle  is  adequate. 

Several  recommendations  for  future  research  can  be  made  concerning 
the  drag  minimization  problem.  The  two  classes  of  bodies  considered  in 
this  study  are  constrained  to  be  well  behaved  according  to  previous 
hydrodynamic  experience.  The  analyses  for  the  various  constraint 
boundaries  are  complicated  and  lack  generality.  It  may  prove  beneficial 
to  develop  a  more  general  class  of  bodies  using  orthogonal  polynomials, 
for  example.  Rather  than  deriving  constraint  boundaries,  it  would  be 
simpler  to  make  a  direct  check  of  the  profile  and  its  derivatives;  how¬ 
ever,  knowledge  of  the  constraint  boundaries  does  provide  insight  into 
the  nature  of  the  problem. 
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The  one  outstanding  weakness  of  the  modified  Complex  Method  is  its 
stopping  condition.  The  problem  is  that  a  large  number  of  function 
evaluations  are  wasted  in  the  process  of  deciding  that  the  best  design 
is  sufficiently  near  a  local  minimum.  A  more  sophisticated  terminating 
strategy  would  be  beneficial. 

In  the  area  of  drag  prediction,  the  major  weakness  appears  to  be 
the  modeling  of  the  transition  region.  The  assumption  of  point  transi¬ 
tion  may  be  poor  at  low  Reynolds  numbers  for  which  there  is  no  abrupt 
increase  in  skin  friction  as  predicted  by  the  drag  model  used  in  this 
study.  The  use  of  the  planar  flow  Michel -e9  correlation  to  predict 
transition  of  ax i symmetric  boundary  layers  probably  introduces  some 
error.  Experimental  verification  of  the  results  of  this  study  is 
certainly  desirable. 

As  far  as  hydrodynamic  design  is  concerned,  the  next  logical  step 
is  to  include  the  propeller  or  propulsive  jet  effects  in  the  drag  model. 
This  would  change  the  problem  from  drag  to  power  minimization.  Since 
the  presence  of  the  propeller  or  propulsive  jet  will  change  the  flow 
field,  it  is  expected  that  the  minimum  power  body  shapes  will  differ 
from  those  for  minimum  drag.  As  mentioned  above,  it  seems  desirable  to 
further  constrain  the  problem  to  avoid  high  sensitivity  to  early  transi¬ 
tion. 
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APPENDIX  A 


DERIVATION  OF  DISPLACEMENT  THICKNESS  FOR  MASS 
CONSERVATION  IN  AN  EXTERNAL  AX  SYMMETRIC  FLOW 


This  appendix  presents  the  derivation  of  the  displacement  thick¬ 
ness  6gX  for  mass  conservation  in  an  external  axisymmetric  flow.  The 
principal  notation  is  shown  in  Figure  Al.  The  zero  incidence  flow  with 


Figure  Al.  External  Axisymmetric  Flow  Notation. 


reference  velocity  is  parallel  to  the  centerline  X-axis.  The  bound¬ 
ary  layer  is  specified  in  terms  of  the  curvilinear  x-y  coordinate 
system.  At  a  point  in  the  boundary  layer  the  velocity  is  u  =  u(x,  y). 
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which  approaches  the  external  value  ue(x)  as  y  approaches  the  boundary 
layer  thickness  6.  The  local  tangent  has  angle  a;  the  local  wall  radius 
is  r0.  It  is  also  convenient  to  use  the  relationship  r  =  r0  +  y  cos  a, 
where  r  is  the  radial  coordinate  of  a  point  in  the  flow.  The  local 
density  p(x,  y)  and  external  density  pe(x )  are  retained  in  the  deriva¬ 
tion  to  follow. 

The  displacement  thickness  concept  equates  the  retarded  boundary 
layer  mass  flow  to  a  displaced  inviscid  mass  flow  of  constant  local 
velocity  ue(x).  This  may  be  written  as 

(Viscous  Mass  Flux)  =  (Inviscid  Mass  Flux) 
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pu  2irr  dy 
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We  assume  that  outside  the  boundary  layer,  i.e.,  for  y  >  6,  that 
pu  =  pe  ue  so  that  equation  (A.l)  reduces  to 
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This  may  be  written  as 
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(A. 4) 


Replacing  r  by  r0  +  y  cos  a  in  the  left-hand  side  of  equation  (A. 4)  and 
dividing  through  by  rc  gives 
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for  r0  /  0.  If  y  cos  a/rQ  «  1 »  then  we  obtain 


Pe  ue 


(A. 6) 


which  is  the  definition  of  s*  as  given  by  equation  (2.17)  where  p  =  pe 
is  constant. 

When  y  cos  a/rQ  cannot  be  neglected,  then  integration  of  the  left 
hand  side  of  equation  (A.5)  gives 


r*  .  cos  a,t£*  \  2  _  f+ 
°ax  +  -g— 1 l®ax>  ’  6 


(A.7) 


where  6*,  as  defined  by  equation  (2.17),  has  replaced  the  integral  on 
the  right-hand  side  of  equation  (A.5).  Equation  (A.7)  is  quadratic  in 
6ax;  applying  the  quadratic  formula  gives 


ax  cos 


M- 


1  +  (l  + 


2  cos 


2-2  S*F ' 
o  J 


(A.8) 


which  is  the  desired  result. 
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APPENDIX  B 

DERIVATION  OF  YOUNG* S  FORMULA  FOR  DRAG  OF 
AXISYMMETRIC  BODIES  AT  ZERO  INCIDENCE 


Young's  formula  is  based  on  the  boundary  layer  momentum  integral 
equation.  If  the  pressure  gradient  across  the  boundary  layer  and  wake 
is  negligible,  the  momentum  integral  equation  for  axisymmetric  bodies 
may  be  written 


ue 


(H  +  2)x 


2irr 


(B.l) 


where  x  is  the  streamwise  coordinate  along  the  body  surface,  ue  is  the 
velocity  at  the  edge  of  the  boundary  layer,  u^  is  du6/dx,  tw  is  the 
shear  stress  at  the  wall,  p  is  the  fluid  density,  r  is  the  radial  co¬ 
ordinate  measured  from  the  body  axis  and  r  =  r0  +  y  cos  a,  where  r0  is 
the  wall  radius,  a  is  the  angle  between  the  surface  tangent  and  the  body 
axis  in  a  meridional  plane,  y  is  the  coordinate  normal  to  the  wall,  x  is 
the  momentum  area,  and  H  is  the  shape  factor. 

For  the  non-separatinq  wake,  the  momentum  integral  equation  is 
assumed  to  apply  and  the  skin  friction  tw  is  zero  so  that  equation  (B.l) 
may  be  written 

&  *  <H  *  2>X  *  0  (B.2) 

where  for  the  wake  region  x  is  the  streamwise  coordinate  measured  from 
the  tail  along  the  axis  of  symmetry,  9y/9x  is  replaced  by  dy/dx,  and  the 
radial  coordinate  r  reduces  to  y. 


The  momentum  area  of  the  wake  is  defined  as 


x  * Zv  r*  °  ■  £)y  dy 


(B.3) 


and  the  displacement  area  of  the  wake  is  defined  as 


A  =  2tt 


C(,-£)ydy 


(B.4) 


so  that  the  shape  factor  is  defined  as 


H  =  A/x 


(B.5) 


Rearranging  equation  (B.2)  and  using  the  fact  that 


d(*n  ue)  = 


(B.6) 


so  that 


±  ( on  He  )  -  d  ,  ue  >  1 

dx  }  dx  (  C  }  uiTDT 


(B.7) 


we  may  write 


1  dx.  _ 

X  dx 


-<H.2)£*n£ 


(B.8) 
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or 


ue 


d(An  x)  =  -  (H  +  2)  d(Jtn  p-  ) 


(B.9) 


Integrating  from  the  tail  downstream  to  infinity  we  obtain 


fOc 

d  An 

J()i 


x 


(H  +  2)  d(An  ^  ) 
'■'00 


(B.10) 


In  X*.  -  in  Xi  =  - 


(H  +  2)  An 

''OO 


0. 


0, 


II 


An 


C 


dH 


(B.l I ) 


where  subscripts  ()w  and  ()l  denote  quantities  at  infinity  downstream 
and  at  the  tail,  respectively,  and  we  have  used  integration  by  parts  on 
the  right-hand  side  of  equation  (B.10).  The  value  of  H  at  ir^inity  is 


unity  [12].  Since 


=  1 ,  we  have  that 


*n  Xoo  = 


(B.12) 


which  can  be  written 


*  X« 

I 

\ 


^e  JH,+2) 


(B.13) 


1-.^ 


7' 


The  exponential  term  may  be  written 


K 


exp  An  77-  dH 

*  1  Uno 


=  exp 


l"1  l"  k 


(B.14) 


At  this  point  Young  assumes  that  An  p  varies  linearly  with  H  so  that 

^00 

the  integral  in  equation  (B.14)  is  approximately  the  area  of  a  triangle. 
This  gives 


fHi  u 

I,  m{€) 


(B.15) 


rH*  u 

U 


, Hrl, 

- 

ue  i 


(B.16) 


Using  equation  (B.16)  in  equation  (B.13)  and  replacing  =  with  =  gives 


x-  *  x.<  Sr  ),  2 


(B.17) 


Equation  (B.17)  relates  the  momentum  area  of  the  wake  at  infinity  x  to 

oo 

the  momentum  area  of  the  wake  at  the  tail  Xi  a°d  other  trailing  edge 


which  is  identical  to  equation  (2.15).  We  have  tacitly  assumed  here 
also  that  the  wake  momentum  area  at  the  tail  is  equal  to  that  of  the 
boundary  layer  there.  This  seems  reasonable  since  the  momentum  defect 
of  viscous  flows  must  change  in  a  continuous  manner. 
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APPENDIX  C 

DETAILED  STRATEGY  OF  MODIFIED  COMPLEX  METHOD 

The  word  flow  chart  given  below  gives  the  detailed  strategy  of  the 
Complex  Method  [22]  as  used  in  the  present  study.  The  strategy  in¬ 
cludes  modifications  cue  to  Guin  [23]  to  cope  with  nonconvex  boundaries 
as  well  as  the  generation  of  new  search  directions  when  the  primary 
direction  fails.  The  starting  procedure  is  a  modification  of  Box's 
original  method. 

A.  Input 

1.  Number  of  independent  parameters  N. 

2.  Number  of  vertii  es  in  complex  figure  K. 

3.  Maximum  number  or  performance  function  evaluations  IPMAX. 

4.  Tolerances  zx  ,  e2,  t3. 

5.  Expansion  factor  a. 

6.  Contraction  factor  B. 

7.  Lower  constant  boundaries  ai_.,  i  =  1,  ...»  N,  used  during 
initial  complex  generation. 

8.  Upper  constant  boundaries  ay. ,  i  =  1 ,  . . . ,  N,  used  during 
initial  complex  generation. 

B.  Initial  Complex  Generation 

1 .  Set  j=l. 

2.  Does  j  exceed  K? 

a.  Yes:  Go  to  B.4. 

b.  No:  Randomly  generate  jth  vertex  using  next  N  elements  of 

'  _  random  number  sequence  indicated  by  (rn-j ) ,  i  =1,  ...,  N. 

f 

r 

f. 

t 


I 

t 


a-j]j  =  (rni)(aui  -  aLi)  ,  0  <  rni  <  1 
i  =  1,  ...,  N 

3.  Is  random  vertex  f[j  feasible? 

a.  Yes:  Evaluate  performance  function  PF-:  for  vertex  a-?. 

Set  j  =  j  +  1 .  J  ”J 

Go  to  B.2. 

b.  No:  Has  random  generation  of  vertex  a*  failed  more  than 

1000  times? 

1)  Yes:  Abort  program. 

CALL  EXIT. 

2)  No:  Try  another  random  generation  of  vertex  a-;. 

Go  to  B.2.b. 

4.  Set  abest  =  best  vertex  of  initial  complex. 

Set  PPbest  =  PF  value  of  2-besf 
Go  to  C.4. 


Search  Procedure 

1 .  Is  newest  vertex  better  than  <LDest? 

a.  Yes:  An  improvement  has  occurred. 

Increment  improvement  counter  IMPRV  =  IMPRV  +  1 
Set  a_best  =  new  vertex. 

Set  PFbest  =  PF  value  for  new  -best' 

Go  to  C. 2. 

b.  No:  Continue. 

2.  Check  stopping  condition:  Are  five  best  PF  values  within 
relative  e3-neighborhood  of 

a.  Yes:  Tentative  convergence.  Are  corresponding  five  best 

vertices  within  relative  e3-neighborhood  of  iLbest? 

1)  Yes:  Convergence  achieved.  Go  to  C. 11. 

2)  Nc:  Try  one  more  vertex  rejection,  i.e.,  if  the 

previous  vertex  rejection  originated  from 
this  point  in  program,  assume  convergence 
and  go  to  C.ll.  Otherwise,  go  to  C.3. 

b.  No:  No  convergence.  Continue. 

3.  Will  next  performance  function  evaluation  exceed  maximum  number 
IPMAX  allowed? 

a.  Yes:  Output  results  obtained  so  far. 

CALL  EXIT. 


IjlUlWAMWW1,  SWBJ 


I 


135 


b.  No:  Continue. 

4.  Reject  worst  vertex.  Rejected  vertex  is  ar. 

5.  Compute  centroid  cgr  of  remaining  vertices. 

1  K 

c9ri  =  ’(k-'IJ  .^-i  ai^j  *  1 •••»  N 
3”  * 
j^r 

6.  Is  cg£  in  feasible  region? 

a.  Yes:  Set  NCGVIO  =  0 

Go  to  C.7. 

b.  No:  Set  NCGVIO  «  1 

7.  Generate  trial  vertex  a_f 
a  t  =  cgr  +  a(cgr.  -  ar) 

8.  Is  trial  vertex  a^  in  feasible  region? 

a.  Yes:  Evaluate  trial  vertex  performance  function  PFt. 

Go  to  C.10. 

b.  No:  Does  NCGVIO  =  1? 

1)  Yes:  Go  to  C.9. 

2)  No:  Is  trial  vertex  a^  within  relative  e2-neigh- 

borhood  of  cgr? 

a)  Yes:  Set  =  cgr  which  is  feasible. 

Evaluate  PF^.  Go  to  C.10. 

b)  No:  Move  trial  vertex  ^  amount  8  to- 

ward  cgr.  New  a^  =  3(old  a_t  +  cgr). 
Go  to  C.8. 

9.  Both  centroid  cgr.  and  trial  vertex  a^  are  not  feasible.  Reject 
entire  complex  figure.  Have  more  than  10  complex  rejections 
occurred? 

a.  Yes:  Abort  search. 

CALL  EXIT. 

b.  No:  Reset  lower  bounds  aL_.  and  upper  bounds  a^j  1  =  1 , 

...,  N,  to  coincide  with  unfeasible  centroid  cgr  and 
best  previous  vertex  abest.  Swap  components  if 

necessary  to  insure  that  ai^  <.  a^*  i  =  1*  ...»  N. 

Go  to  B.l . 

10.  Is  trial  vertex  performance  function  PFt  better  than  second 
worst  PF  value? 


i 


iflagafririitefl 
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a.  Yes:  Replace  worse  ve^ex  with  vertex  a^. 

Go  to  C. 1 . 

b.  No:  Is  trial  vertex  within  relative  e2-reighborhood  of 

oentroid  egr? 

1)  Yes:  Try  a  new  search  direction.  Reject  the  next 

worst  vertex  and  retain  previously  rejected 
vertex.  Rejected  vertex  is  a_r.  Have  all  K 
vertices  been  successively  rejected? 

a)  Yes:  Abort  search. 

CALL  EXIT. 

b)  No:  Go  to  C.5. 

2)  No:  Move  trial  vertex  a^  amount  6  toward  centroid 

cqr.  New  a_t  =  8(old  a^.  +  '.^r).  Go  to  C.8. 

11.  Output  optimum  parameters  a*. 

‘  -best 


12.  END 


APPENDIX  D 


DETAILED  STRATEGY  OF  POWELL'S  METHOD  OF  CONJUGATE  DIRECTIONS 


The  word  flow  chart  given  below  presents  the  strategy  for  Powell's 
Method  of  Conjugate  Directions  [25].  The  only  modification  from  Pow- 
ell's  original  procedure  is  the  procedure  for  computing  the  new  step 
size  STEP  used  in  the  linear  search  routine;  see  item  B.15  in  the  out¬ 
line  below.  Following  Powell's  method  is  a  word  flow  chart  for  the 
linear  search  strategy  involving  a  parabolic  interpolation  scheme. 


D.l  Strategy  for  Powell's  Method 

A.  Input 

1.  Number  of  independent  parameters  N. 

2.  Feasible  initial  guess  vector  a0  =  (a^,  a^,  ...»  a0^). 

3.  Set  of  N  linearly  independent  normalized  search  direction 
vectors  (£_j»  fj_2>  •••>  JiN^* 

4.  Lower  and  upper  scaling  vectors,  ay  and  a^y  ,  respectively, 
where  a_y  -  (aj_^,  ay^»  •••*  ay^)  and  a_y  -  (ay^>  ay^,  ...,  ay^). 

5.  Maximum  number  of  search  cycles  ICYMAX. 

6.  Initial  scaled  step  size  STEP  used  in  linear  search  routine. 

7.  Convergence  tolerance  e3. 

B.  Search  Procedure 

1.  Initialization. 

a.  Scale  initial  guess  vector  to  X^0  ,  where  X0i  =  (a0i  -  a^) 
/(ay-  -  a l ^ ) »  i  _  1*  •••»  N. 


*ui 


ift'iYkw 


b.  Evaluate  performance  function  PF0  of  initial  guess. 

2.  Start  search  cycle. 

a.  Set  maximum  PF  change  Am  =  0. 

b.  Sec  direction  number  For  maximum  PF  change  im  =  1. 

c.  Set  direction  index  j  =  1. 

3.  CALL  COGGIN  (linear  search  routine):  from  base  point  X.-_i 

J  • 

perform  linear  minimization  along  direction  vector  Find 
the  minimum  point  Xj  and  its  performance  function  value  PFj. 

4.  Test  for  maximum  PF  change.  Is  |PFj  -  PFj_-j  |  greater  than 


Am 

? 

a. 

Yes: 

Set  4m  -  |PFj  -  PFj.,1. 
Set  im  =  j.  Go  to  B.5. 

b. 

No: 

Go  to  B.5. 

Increment  direction  index  j  =  j  +  1. 

a. 

Yes: 

Continue  cycle,  go  to  B.3. 

b. 

No: 

Cycle  is  complete. 

Go  to  B.6. 

6.  Set  Xbest  XN  and  PPbest  =  pFN* 

7.  Compute  new  direction  vector  ^  =  (Xj»j  -  \0)t  that  is,  y-j  = 

XNi  ~  Xq.j  i  1  1,  « « « ,  N. 

8.  Compute  trial  point  =  Xpj  +  that  is,  y^  =  2Xn^  -  X0^, 

i  =  1,  ...,  N.  Is  -  X0|  less  than  ( .001 ) (e3 )? 

a.  Yes:  Convergence  likely.  Go  to  B.14. 

b.  No:  Compute  PFt  =  PF(^t).  Go  to  B.9. 

9.  Perform  Powell's  first  inequality  check  on  new  direction.  Is 
(PFt  -  PF 0)  ^  0? 

a.  Yes:  New  direction  is  not  promising. 

Go  to  B.13. 

b.  No:  Go  to  B.10. 

10.  Perform  Powell's  second  inequality  check  on  new  direction.  Is 
(PF0  -  2PF„  +  PFt)(PF0  -  PFN  -  im)  >  i-  An(PF0  -  PFt)? 

a.  Yes:  New  direction  is  not  promising. 

Go  to  B.13. 

b.  No:  Go  to  B.ll . 
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11.  New  direction  is  promising.  Remove  direction  vector  of  maxi¬ 
mum  PF  change,  .  Put  new  direction  vector  in 

a.  For  i  =  1,  ...»  im  -  1 :  new  =  old  £j. 

b.  For  i  =  im  -  1 ,  . .  • ,  N  -  1 :  new  £.j  =  old  . 

c.  For  i  =  N:  new  =  h/Ih.1* 

12.  CALL  COGGIN  (linear  search  routine):  from  base  point 
perform  linear  minimization  along  new  direction  vector 
Find  new  minimum  point  t  and  its  performance  function 

value  PFj)est.  Go  to  B.14. 

13.  New  direction  is  not  promising.  Retain  old  search  directions: 

a.  For  i  =  1 ,  ...»  N:  new  =  old  £.j . 

b.  Go  to  B.14. 

14.  Test  for  convergence, 
a.  Is  PF0  -  PFbe$t 

P'P0 .  s' 

1)  Yes:  Go  to  B.14.b. 

2)  No:  Go  to  B.15. 


b.  Is 


X°i  "  Xbest  i 


*  E3? 


i  =  1 ,  . . . ,  N 

1)  Yes:  Convergence  achieved. 

Go  to  B.17. 

2)  No:  Go  to  B.15. 

15.  No  convergence  yet.  Compute  new  step  size  STEP: 


a. 

AX  = 

|  -best  "  -o 

b. 

apf  = 

|PFbest  -  PF°| 

c. 

STEP 

=  /SPF. 

d. 

Is  ST 

EP  <  (0.1 ) ( AX ) ? 

1)  Yes:  Set  STEP  =  (0.1)(AX). 

Go  to  B.15.e. 

2)  No:  Go  to  B.15.e. 


■v  **  «•%.  'fiys 


e.  Is  STEP  >  (0 . 5) (old  STEP)? 

1)  Yes:  Set  STEP  -  (0.5)(old  STEP). 

Go  to  B.15.f . 

2)  No:  Go  to  B.15.f. 

f-  Se‘io=Ibest.PFo  =  PFbest. 

16.  Will  next  cycle  exceed  maximum  number  ICYMAX  of  cycles  allowed? 

a.  Yes:  Output  results  obtained  so  far. 

CALL  EXIT. 

b.  No:  Go  to  B.2. 

17.  Convergence  achieved. 

a.  Scale  optimum  X*  back  to  a*:  a*  =  a^.  +  X*  (a^  -  a^), 
i  =  1,  ....  N. 

b.  CALL  EXIT. 

18.  END. 


Strategy  for  Linear  Search  Routine 


The  word  flow  chart  below  gives  the  detailed  strategy  for  the 
parabolic  interpolation  scheme  used  for  the  linear  search  in  the  pres¬ 
ent  study.  The  method  always  brackets  the  minimum  along  the  line  of 
search  before  applying  the  parabolic  interpolation. 


A.  Input  from  Calling  Program  POWELL 

1.  Current  base  point  is  £  with  performance  function  value  FX. 

2.  Current  step  vector  is  S_  =  STEP  £,  where  £  is  the  current 

normalized  search  direction  and  S  =  (St ,  S2,  ....  Sfj).  S5  = 
(STEP) (Si ) .  i  =  1,  ...,  N. 

B.  Linear  Search  Procedure. 

1.  Initialization. 

a.  Set  step  multiplier  D  =  1. 

b.  Set  distance  DA  =  0  for  point  A. 

c.  Set  distance  DB  =  0  for  point  B. 
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d.  Set  distance  DC  =  0  for  point  C. 

e.  Set  performance  function  FA  =  FX  fo>*  point  A. 

f.  Set  performance  function  FB  =  0  for  point  B. 

g.  Set  performance  function  FC  =  0  for  point  C. 

h.  Set  step  counter  K  =  -2. 

i.  Set  linear  search  convergence  tolerance  TOL  =  3e3. 

2.  Start  linear  search. 

a.  Compute  point  y.  =  X.  +  D!S,  where  y-j  =  Xj  +  D  S j ,  i  =  1,  ...»  N. 

b.  Compute  performance  function  F  for 

c.  Increment  step  counter  K  =  K  +  1. 

3.  Is  F  >  FA? 

a.  Yes:  A  bracket  point  has  been  found. 

Go  to  B.4. 

b.  No:  Performance  function  is  still  decreasing. 

1)  Reset  points  A,  B,  and  C 

DC  =  DB  FC  =  FB 

DB  =  DA  FB  =  FA 

DA  =  D  FA  =  F 

2)  Compute  new  step 

For  D  >  0:  new  D  =  2(old  D)  +  1 . 

For  D  <  0:  new  D  =  2(old  D)  -  1 . 

3)  Go  to  B.2. 

4.  A  bracket  point  has  been  found.  Is  K  >.  0? 

a.  Yes:  Both  bracket  points  have  been  found. 

Set  up  points  A,  B,  and  C  so  that  the  minimum  point 
B  is  bracketed  by  the  points  A  and  C: 

DC  =  DB  FC  =  FB 

DB  =  DA  FB  =  FA 

DA  =  D  FA  =  F 

Go  to  B.5. 

b.  No:  First  bracket  point  has  been  found  after  first  step 

along  line  of  search.  Reverse  direction  and  continue. 

1)  Set  print  B 

FB  =  F  DB  =  D 

2)  Change  sign  of  step  multiplier 
new  D  =  -  old  D 

3)  Go  to  B.2. 


5.  Proceed  with  parabolic  interpolation. 

a.  Compute  location  D  of  minimum  on  parabolic  arc. 

b.  Is  point  D  between  DA  and  DC? 

1)  Yes:  Go  to  B.5.c. 

2)  No:  Parabolic  interpolation  has  failed.  Use  best 

point  B  as  local  minimum. 

Go  to  B.6. 

c.  Perform  parameter  convergence  check: 

Is  |D  -  DA |  /  | DA  -  DC {  <  TOL? 

1)  Yes:  Convergence  achieved.  Use  point  B  as  local 

minimum.  Go  to  B.6- 

2)  No:  Go  to  B.5.d. 

d.  Perform  function  convergence  check.  Compute  point  =  X.  + 
D  and  its  performance  function  value  F.  Is 

| FB  -  F|  /  | FB |  <  TOL? 

1)  Yes:  Convergence  achieved.  Use  smaller  of  FB  or  F  as 

local  optimum. 

Go  to  B.6. 

2)  No:  No  convergence  yet.  Reset  points  so  that  minimum 

point  B  is  bracketed  by  points  A  and  C. 

Go  to  B.5. 

6.  Convergence  to  local  minimum  along  line  of  search  has  been 
achieved.  Compute  approximate  second  derivative  associated 
with  present  search  direction,  32PF/3£2. 

7.  RETURN  TO  POWELL. 
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APPENDIX  E 

DERIVATION  OF  PARAMETRIC  BODY  PROFILES 

This  appendix  presents  the  derivations  for  the  parametrically  de¬ 
fined  body  profiles  considered  in  the  present  study.  The  pertinent 
constraint  boundaries  are  also  derived.  The  procedures  follow  those 
reported  by  Granville  (1969)  [29]  in  which  the  body  is  divided  into 
sections  at  convenient  axial  locations,  e.g.,  maximum  diameter  point. 

E.l  Rounded-Nose  Forebody  Section  [29] 

The  rounded-nose  forebody  and  its  dimensional  pe-^meters  are  showw 
in  Figure  El.  The  maximum  diameter  D  occurs  at  the  axial  location 
the  curvature  at  this  point  is  .  The  nose  radius  of  curvature  is  Rn. 


Figure  El.  Rounded-Nose  Forebody  Section. 
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The  boundary  conditions  for  this  profile  Y(X)  are  listed  below  where 
primes  on  Y(X)  denote  differentiation  with  respect  to  X: 


1.  Y(0)  =  0 


(E. la) 


2.  Y'(0)  =  »  or  dX/dY]„3Q  =  0 


(E. lb) 


3.  l/(d2X/dY2)]x=Q  =  Rn 


4.  Y(Xm)  =  D/2 


5.  Y' (Xm)  =  0 


(F.lc) 


(E.ld) 


(E.le) 


6.  Y"(Xm)  =  Ki 


(E.lf) 


where  Rn,  D,  Ki,  and  Xm  may  all  vary. 

It  is  possible  to  reduce  Y(X)  to  a  nondimensional  profile  with 
only  two  variable  boundary  conditions.  We  define  the  nondimensional 
profile  to  be 


x  =  X/XM 


( E . 2a ) 


y  =  2Y/D 


(E.2b) 


so  that  the  boundary  conditions  given  by  equations  (E.l)  become 


1.  y(0)  =  0 


(E.3a) 


2.  dx/dy]x_Q  =  0 


(E.3b) 


3.  l/(d2x/dy2)]x=0  =  4Xm  Rn/D2  = 


(E.3c) 


4.  y(l)  =  1 


(E.3d) 


5.  y'O)  =  0 


/'’■'7V  JTf? »T'*,rv7T?’T»'’KY> I"* 


(E.3e) 


6.  y"(l)  =  2X*  Kj/D  =  -k, 


(E.3f) 


where  primes  on  y(x)  denote  differentiation  with  respect  to  x.  The  var¬ 
iable  boundary  conditions  are  rn  and  kx ;  two  free  parameters  make  it 
possible  to  conveniently  plot  constraint  boundaries  for  the  forebody. 

The  intent  here  is  to  derive  the  simplest  polynomial  expressions 
which  satisfy  the  boundary  conditions  given  by  equations  (E.3).  Fur¬ 
thermore,  we  shall  attempt  to  make  the  expressions  linear  in  the  free 
parameters  rn  and  kx.  For  this  purpose  we  define  the  "quadratic"  poly¬ 
nomial 


f(x)  =  y2(x) 


n 

y  bi  xi 

i-O 


(E.4) 


where  the  bj,  i  =0,  ....  n,  are  to  be  determined  using  the  boundary 
conditions  given  by  equations  (E.3).  Differentiating  equation  (E.4) 
with  respect  to  x  gives 


f'  =  ^-y2  =  2yy'  =  J  (i)  ^  x1’1 


(E.5) 


f"  =  TT  y2  =  2yy"  +  2y'2  =  l  (i)(i-l)  bj  x1_?- 
dx*  i=? 


(E.6) 


f,M  =  — r  y2  =  2yy"'  +  6y"y’ 
dx3 


l  (i)(i-D(i-2)  bi  xi_3  (E.7) 
i=3 
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where  primes  on  f(x)  denote  differentiation  with  respect  to  x.  Equa¬ 
tion  ( E . 5)  may  be  written  as 

2y  '  (  Ji(1)  b)  x1'i  J  37  (E’8) 

and  differentiating  equation  (E.8)  with  respect  to  y  gives 


2 


l  (0(1-1 ) 


1=2 


(E.9) 


For  the  case  x  =  0,  equations  (E.4),  (E.8),  and  (E.9)  reduce  respective¬ 
ly  to 


y(o)  =  o  =  b0 


(E.10) 


2y(0)  «  0  =  b,  [dx/dy]xa() 


(E.ll) 


and 


2  =  2b2  [dx/dy]’_0  +  b,  [d2x/dy2]x=0 


(E.12) 


where  equation  (E.10)  has  been  inserted  into  equation  (E.11).  Equation 
(E.l 1 )  implies  that  bt  or  [dx/dy]x_Q  or  both  equal  zero.  If  bx  equals 

zero,  then  boundary  condition  (E.3c)  cannot  be  satisfied  by  equation 
(E.12).  Hence  bj  f  0  and  boundary  condition  (E.3b)  is  automatically 
satisfied  when  y (0)  =  0.  This  is  the  motivation  for  selecting  the 
"quadratic"  polynomial,  equation  (E.4),  at  the  outset.  From  equations 
(E.10),  (E.ll ) ,  and  (E.12)  we  obtain  Immediately  that 


[  Using  equations  ( E . 1 0)  and  (E.13),  the  boundary  conditions  (E.3d,  e5  f) 

\  give 


2rn  +  .1  bi  =  1 
i=2 


2rn  +  HDbi  =  0 
i=2 


( E. 14a) 


(E . 14b) 


l 

l  ?  (i)(i-l)bi  =  -2k,  ( E. 14c) 

r  i=2 

l 

j 


Setting  -  =  4  yields  a  unique  solution  for  equations  (E.14). 

f  Rather  than  solving  for  the  b,,  i  =  2,  3,  4,  by  direct  elimination, 

; 

I'  we  follow  Granville  [29]  and  postulate  the  existence  of  a  function  y2(x) 

which  is  linear  in  rn  and  kx.  This  is  justified  since  y2  is  linear  in 
the  bj,  i  =  2,  3,  4,  which  are  linear  in  rn  and  k,.  Thus  we  write 


1 


y2(x)  =  rn  Fjtx)  +  kjF2(x)  +G(x) 


(E.15) 


c 

t 


i 


£ 

! 


if 

sf 

I 


where  F^x),  F2(x),  and  G(x)  are  polynomials  of  degree  n  =  4,  and  rn  and 
k,  are  independent  and  arbitrary  for  the  moment.  The  six  boundary  con¬ 
ditions  (E.3)  imply  the  following  five  boundary  conditions  for  F,(x), 
F2(x),  and  G(x)  for  arbitrary  values  of  rn  and  k,: 

1.  y (0)  =  0  implies  that  F,(0)  =  F2(0)  =  G(0)  =  0 

2.  y{l)  =  1  implies  that  F,(l)  =  F2(l)  =  0  and  G(l)  =  1 


■J- - - - 
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3.  dx/dyJx=0  =  0  implies  that 

2y(0)y'(0)  =  bx  =  2rn  =  rnF;  (0)  +  k&  (0)  +  6'  (0) 

so  that  Fj  (0)  -  2  and  F2'(0)  =  G'(0)  =  0 

4.  y'(l)  =  0  implies  that  Fj'(l)  =  F2(l)  =  G ' (1 )  =  0 

5.  y"(l)  =  -kj  implies  that 

-2k!  =  rnF" (L)  +  kjF2"(l )  +  G"(l)  so  that 

F” ( 1 )  =  G"(l )  =  0  and  F2"(l)  =  -2 

Knowing  the  polynomial  degree  and  the  boundary  conditions,  one  can 
v/rite  down  the  functions  Fj(x),  F2(x),  and  G(x)  almost  by  inspection. 

For  Fj(x)  we  have  that 

F,(0)  =  0  FjO)  =  0 

Fx{0)  =  2  Fj(l)  =  0 

FKD  =  0 

so  that  we  may  write  immediately  that 
F^x)  =  Cj  x(x  -  l)3 

to  satisfy  the  homogeneous  boundary  conditions.  Applying  the  final 
condition  FJ{0)  =  2  gives 

Fx(x)  =  -2x(x  -  l)3  (E.16) 

'  .  i lari y  we  may  write 

F2(x)  =  -x2(x-l)2  (E.17) 

For  G(x)  we  h?*  e  that 


*.  V. 
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G(0)  =  0  G(1 )  =  1 

G'(0)  =  0  G ' (1 )  =  0 

G"(l )  =  0 


so  we  may  write 


^  G(x)  =  Cj  x(x  -  l)2 


and  integrating  gives 

G(x)  =  c,(ix4  -  fx3  +|x2)  +  c2 


For  G(0)  =  0  we  have  c2  =  0  and  for  G(l)  =1  we  obtain  finally  that 

G(x)  =  x2(3x2  -  8x  +  6)  (E.18) 


We  now  have  completely  defined  y2{x)  by  equations  (E.15),  (E.16),  (E.17), 
and  (E.18)  which  satisfy  the  boundary  conditions  (E.3).  Inserting  these 
results  into  equation  (E.2b),  and  using  the  fact  that  fr  =  L/D,  we  ob¬ 
tain  the  forebody  profile  equations  (4.3)  and  (4.4)  in  Chapter  4. 

Constraint  Boundaries  for  rn  and  kj.  From  physical  considerations 
we  require  a  non-negative  radius  of  curvature  Rn  at  the  nose  and  no.:~ 
positive  curvature  Ka  at  the  maximum  diameter  section.  We  have  defined 
rn  and  kx  so  that  rn  >  0  and  kj  >  0  satisfy  these  physical  requirements. 

Mainly  because  of  previous  hydrodynamic  experience,  an  additional 
requirement  is  introduced,  namely,  that  of  no  inflection  point  on  the 
forebody  [29].  The  limiting  case  occurs  when  the  second  derivative 
touches  zero  but  does  not  change  sign  and  is  expressed  as 


d2Y(X)/dX2  =  0 


(c .19a) 


^ .MU. 
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d3Y(X)/dX3  =  0  (E .! 9b) 

for  0  ^  X  i  Xm.  If  a  solution  exists  for  equations  (E.19),  it  repre¬ 
sents  the  rn  versus  kx  curve  for  which  a  limiting  inflection  point 
occurs  somewhere  on  0  <  X  <  Xm. 

Using  equations  (E.2),  (E.5),  (E.6),  (E.7),  and  (E.15),  equations 
(E.19)  reduce  to 

2ff"  -(f)2  =  0  ( E . 20a ) 

=  0  (E.20b) 

The  procedure  is  to  fix  the  location  x  on  0  <  x  <  1  where  the  limiting 
inflection  point  occurs  and  to  solve  equations  (E.20)  for  the  corre- 
ponding  values  of  rn  and  kx.  In  this  s^nse  rn  and  are  related 
through  the  parameter  x.  The  Secant  Method  [30]  is  used  to  obtain  the 
solution  which  is  given  in  Table  El  and  plotted  in  Figure  20  in  Chap¬ 
ter  4. 

E.2  Pointed  Aftbody  Section  [29] 

The  pointed  aftbody  and  its  dimensional  parameters  are  shown  in 
Figure  E2.  The  maximum  diameter  D  occurs  at  the  axial  location  Xm; 


i,  Y . 


Figure  E2.  Pointed  Aftbody  Section. 
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Table  El.  Solution  for  Limiting  Inflection  on  Forebody. 


X 

*1 

rn 

X 

ki 

rn 

.000 

4.0000 

.0000 

.520 

3.7301 

1.1187 

.020 

4.0398 

.0005 

.540 

3.5616 

1.2439 

.040 

4.0794 

.0024 

.560 

3.3711 

1 .3743 

r>cn 

•  vw 

4.1184 

.0057 

.580 

3.1590 

1.5080 

.080 

4.1566 

.0105 

.600 

2.9269 

1.6423 

.100 

4.1940 

.0172 

.620 

2.6774 

1 .7742 

.120 

4.2300 

.0258 

.640 

2.4146 

1.9002 

.140 

4.2646 

.0366 

.660 

2.1436 

2.0166 

.160 

4.2972 

.0499 

.680 

1.8710 

2.1198 

.180 

4.3276 

.0659 

.700 

1.6034 

2.2069 

.200 

4.3552 

.0849 

.720 

1 .3479 

2.2756 

.220 

4.3796 

.1072 

.740 

1.1104 

2.3251 

.240 

4.4000 

.1333 

.760 

.8957 

2.3554 

.260 

4.4160 

.1634 

.780 

.7068 

2.3681 

.280 

4.4266 

.1980 

.800 

.5448 

2.3654 

.300 

4.4311 

.2375 

.820 

.4094 

2.3498 

.320 

4.4285 

.2823 

.340 

.2988 

2.3243 

.340 

4.4178 

.3331 

.860 

.2108 

2.2913 

.360 

4.3977 

.3901 

.880 

.1424 

2.2533 

.380 

4.3670 

.4539 

.900 

.0909 

2.2121 

.400 

4.3243 

.5250 

.920 

.0534 

2.1692 

.420 

4.2681 

.6038 

.940 

.0276 

2.1258 

.440 

4.1968 

.6906 

.960 

.0112 

2.0828 

.460 

4.1088 

.7855 

.980 

.0026 

2.0407 

.480 

4.0027 

.8887 

1.000 

.0000 

2.0000 

.500 

3.8768 

1.0000 
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the  curvature  at  this  point  is  Kr  The  finite  slope  at  the  tail  is  St. 
The  overall  body  length  is  L.  The  boundary  conditions  for  this  profile 
are  listed  below  where  primes  on  Y(X)  denote  differentiation  with  re¬ 
spect  to  X: 


1. 

Y(Xm)  =  D/2 

(E.21a) 

2. 

Y‘(Xm)  =  0 

( E . 21 b) 

3. 

Yn(Xm)  = 

(E . 21c ) 

4. 

Y(L)  =  0 

(E.21d) 

5. 

Y'(L)  =  St 

(E.21e) 

where  D,  Kx,  St,  Xm,  and  L  may  all  vary. 

It  is  possible  to  reduce  Y(X)  to  a  nondimensional  profile  with 
only  two  variable  boundary  conditions.  We  define  the  nondimensional 


profile  to  be 

x  =  (L-X)/(L-Xm)  (E.22a) 

y  =  2Y/D  (E.22b) 

so  that  the  boundary  conditions  given  by  equations  (E.21)  become 

1.  y(0)  =  0  (E.23a) 

2.  y'(0)  =  ~2{L-Xm)St/D  -*  st  (E.23b) 

3.  y(l )  =  1  (E.23c) 

4.  y' (1 )  =  0  (E.23d) 

5.  y"(l )  =  2{L  -  Xm)2  Kj/D  =  -klg  (E.23e) 


where  primes  on  y(x)  denote  differentiation  with  respect  to  x,  From 
equations  (E.3f)  and  (E.23e)  it  is  apparent  that 


^  iC  ( L-Xfn  )/Xjp] 2  =  k1[(l-xm)/xm]‘' 


(E.24) 


where  xm  -  Xm/L.  The  variable  boundary  conditions  are  st  and  klfl;  two 

free  parameters  make  it  possible  to  conveniently  plot  constraint  bound¬ 
aries  for  the  pointed  aftbody. 

We  follow  the  same  procedure  used  in  Section  E.l  for  the  forebcdy 
so  that  equations  (E.4),  (E.5),  (E.6),  and  (E.7)  apply.  For  the  case 
x  =  0,  equations  (E.4)  and  (E.5)  reduce  to 


y(o)  =  0  =  b0 


(E.25) 


2y(0)y'(0)  =  b1 


(E.26) 


In  order  to  satisfy  the  finite  slope  requirement  when  using  the  "quad¬ 
ratic"  polynomial,  equation  (E.4),  then 


bi  =  o 


(E.27) 


so  that  y'(0)  is  indeterminant  in  equation  (E.26).  Setting  bx  =  0  and 
dividing  through  by  2y(x)  in  equation  (E.5),  we  may  apply  L'Hopital's 
Rule  in  the  limit  as  x  approaches  zero  so  that 


lim  y'(x) 
x->0 


l  (i)(i-l)bi  x1’"2 
i=2 


=  st 


(E.28) 


Using  equations  (E.25),  (E.27),  and  (E.28),  the  boundary  conditions 
(E.23c,  d,  e)  give 

2  n 

St  +  I  bi  =  1  (E.29a) 

i=3 

2s|  +  I  (i)b-j  =  0  (E.29b) 

i=3 

2s\  *  j  (1)(i-1)b,  =  -2k,a  ( E . 29c ) 

l  “3 

Setting  n  =  5  yields  a  unique  solution  for  equations  (E.29). 

Rather  than  solving  for  the  b-j,  i  =  3,  4,  5,  by  direct  elimination, 
we  follow  the  procedure  given  in  Section  E.l.  Since  y2(x)  is  linear  in 
the  b-j,  i  =  3,  4,  5,  which  are  linear  in  s\  and  klfl,  we  postulate  that 

y2(x)  =  s\  Fj(x)  +  klg  F2(x)  +  G(x)  (E.30) 

where  F^x),  F2(x),  and  G(x)  are  polynomials  of  degree  n  =  5,  and  s\ 
and  kla  are  independent  and  arbitrary  for  the  moment.  Following  the 

procedures  of  Section  E.l,  we  obtain  finally  the  pointed  aftbody  profile 
equations  (4.5)  and  (4.6)  in  Chapter  4.  In  equation  (4.5),  parameter 
klfl  has  been  eliminated  using  equation  (E.24). 


Constraint  Boundaries  for  s^  and  klg.  From  physical  considerations 

we  require  a  nonpositive  slope  S^  at  the  tail,  and  a  nonpositive  curva¬ 
ture  Kx  at  the  maximum  diameter  section.  We  have  defined  s^  and  klg  so 

that  Sf  >.  0  and  k.  >.  0  satisfy  these  physical  requirements. 

w  *  a 
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Either  no  or  one  inflection  point  on  the  pointed  aftbody  is  al¬ 
lowed.  We  consider  the  noninflected  case  first,  proceeding  as  in 
Section  E.l  for  the  forebody.  The  limiting  inflection  point  condition 
is  expressed  by  equation  (E.19)  for  Xm  <.  X  <  L.  Using  equations  (E.22), 
(E.29),  (E.4),  (E.5),  (E.6),  and  (E.7),  equations  (E.19)  reduce  to  equa¬ 
tions  (E.20).  These  are  solved  using  the  Secant  Method  [30].  The 
solution  is  given  in  Table  E2  and  is  plotted  in  Figure  19  in  Chapter  4. 

Although  Table  E2  contains  an  entry  for  x  =  0,  there  is  actually  a 
singularity  there  with  an  infinite  number  of  solutions.  For  x  =  0, 
equation  (E.20a)  is  identically  satisfied  and  equation  (E.20b)  reduces 
to 


3st  +  kla  =  10 


(E .31 ) 


where  equation  (E.30)  has  been  used.  This  straight  line  is  plotted  in 
Figure  21  in  Chapter  4. 


E.3  Midbody  Section 


The  midbody  section  used  with  the  eight-parameter  tailboom  body  is 
shown  in  Figure  E3  with  its  dimensional  parameters.  The  maximum  diam¬ 
eter  D  occurs  at  the  axial  location  Xm;  the  curvature  at  this  point  is 


Y 


K! 


Figure  E3.  Midbody  Section. 
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Table  E2.  Solution  for  Limiting  Inflection  on  Pointed  Aftbody. 


X 

k  * 

Kla 

st 

X 

k  * 

K*a 

2 

H 

.000 

5.0000 

1.6666 

.520 

7.1756 

3.1789 

.020 

5.0546 

1.6490 

.540 

7.1986 

3.5385 

.040 

5.1118 

1.6317 

.560 

7.1927 

3.9549 

.060 

5.1718 

1.6150 

.580 

7.1518 

4.4330 

.080 

5.2348 

1.5990 

.600 

7.0690 

4.9769 

.100 

5.3008 

1.5840 

.620 

6.9370 

5.5894 

.120 

5.3701 

1.5703 

.640 

6.7480 

6.2706 

.140 

5.4428 

1.5582 

.660 

6.4944 

7.0169 

.160 

5.5191 

1.5482 

.680 

6.1695 

7.8194 

.180 

',.5990 

1.5407 

.700 

5.7690 

8.6614 

.200 

5.6826 

1.5365 

.720 

5.2923 

9.5174 

.220 

5.7701 

1.5362 

.740 

4.7449 

10.3520 

.240 

5.8614 

1.5407 

.760 

4.1403 

11.1212 

.260 

5.9564 

1.5510 

.780 

3.5005 

11.7770 

.280 

6.0549 

1.5685 

.800 

2.8551 

12.2756 

.300 

6.1568 

1.5946 

.820 

2.2371 

12.5868 

.320 

6.2614 

1.6311 

.840 

1.6772 

12.7019 

.340 

6.3682 

1.5799 

.860 

1.1976 

12.6356 

.360 

6.4762 

1./434 

.880 

.8090 

12.4200 

.380 

6.5845 

1.8243 

.900 

.bill 

12.0961 

.400 

6.6914 

1.9257 

.920 

.2954 

11.7039 

.420 

6.7952 

2.0508 

.940 

.1493 

1 » .2769 

.440 

6.8937 

2.2036 

.960 

.0595 

10.8403 

.460 

6.9843 

2.3881 

.980 

.0133 

10.4112 

.480 

7.0640 

2.6089 

1.000 

.0000 

10.0000 

.500 

7.1291 

2.8708 

*  k  = 

Kla 

k, 

*m 

as  given  by  equation 

(E .24) 
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K,.  An  Inflection  point  occurs  at  X^;  at  this  point  the  slope  is  S-j 
and  the  profile  radius  is  R-j.  The  boundary  conditions  for  this  profile 
Y(X)  are  listed  below  where  primes  on  Y(X)  denote  differentiation  with 
respect  to  X: 


1.  Y(Xm)  -  D/2 

(E.32a) 

2.  Y'  (Xjn)  =  0 

(E.32b) 

3.  Y"(Xm)  =  Kx 

(E.32c) 

4.  Y(Xi)  =  Ri 

(E.32d) 

5.  Y'(Xi)=  Si 

(E.32e) 

6.  Y"(Xi)  =  0 

(E.32f) 

where  D,  Kx,  Ri,  Si,  Xm,  and  Xi  may  all  vary. 

It  is  possible  to  reduce  Y(X)  to  a  nondimensional  profile  with  only 
two  variable  boundary  conditions.  We  define  the  nondimensional  profile 

to  be 


x  =  (Xi  -  X)/(Xi  -  Xm)  (E.33a) 

y  =  (Y  -  R1)/(D/2  -  Ri)  (E.33b) 

so  that  the  boundary  conditions  given  by  equations  (E.32)  become 

1.  y (0)  =  0  (E.34a) 

2.  y'(0)  =  -{Xi  -  Xm)Si/(D/2  -  Ri)  =  Si  (E.34b) 

3.  y"{0)  =  0  ( E ,34c) 


- - - -  -  - 
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4.  y(l)  =  1  (E.34d) 

5.  y'(l)  =  0  (E.34e) 

6.  y"(l)  =  (Xi  -  Xm)2K1/(D/2  -  Ri)  =  -klm  (E.34f) 

where  primes  on  y(x)  denote  differentiation  with  respect  to  x.  From 
equations  (E.3f)  and  (E.34f)  it  is  apparent  that 

kim  =  k1(Xi/Xm  -  1)2/(1  -  2Ri/D)  =  kjtxi/Xjj,  -  1)2/<1  -  r-j )  (E.35) 

where  xm  =  Xm/L,  xi  =  X-j/L,  and  r]  =  2R^/D.  The  variable  boundary 
conditions  are  s^  and  klm;  two  free  parameters  make  it  possible  to 

conveniently  plot  constraint  boundaries  for  the  midbody  section. 

Deviating  somewhat  from  Granville's  procedure,  we  postulate  im¬ 
mediately  that 

y(x)  =  klm  F j (x)  +  Si  F2(x)  +  G(x)  (E.36) 


where  Fi(x),  F2(x),  and  G(x)  are  polynomials  of  degree  five  since  there 
are  six  boundary  conditions  given  by  equations  (E.34).  It  is  emphasized 
that  equation  (E.36)  is  an  "ordinary"  polynomial  involving  y(x)  rather 
than  the  "quadratic"  polynomial  used  by  Granville  which  involves  y2(x). 
If  equation  (E.36)  is  an  incorrect  postulate,  then  some  of  the  boundary 
conditions  (E.34)  will  not  be  satisfied. 

For  arbitrary  values  of  klm  and  Si  the  boundary  conditions  for  the 

functions  Fx(x),  F2(x),  and  G(x)  are 


F,(0)  -  0 

F,(l) 

F{(0)  =  0 

F|(l) 

F'i'{0)  =  0 

F'i(l) 

0 

0 

-1 


i 


For  the  function  F:(x)  we  may  write 

Fj(x)  =  c j  x3(x  -  l)2 

which  immediately  satisfies  the  homogeneous  boundary  conditions.  The 
nonhomogeneous  condition  gives 

F,(x)  =  4xJ(x-l)2  (E.37) 

For  the  function  F2(x)  we  may  write 

5 

F2(x)  -  Jbi  x1 

since  the  boundary  conditions  cannot  be  satisfied  by  the  simple  form 
(constant)  xm  (x  -  l)n.  By  direct  elimination  we  obtain 

F2(x)  =  x  -  x3(3xz  -  8x  +  6)  (E.38) 

For  the  function  G(x)  we  may  write 

;j7  G<x)  =  ci  x2(x  ”  ])2 
Integrating  and  applying  G(l)  =1  gives 


'"■/■if*'*' 


nW*V*^*. 
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G(x)  =  x3(6x2  -  15x  +  10)  (E.39) 

Since  all  the  y(x)  boundary  conditions  are  satisfied,  then  apparently 
the  postulated  form  given  by  equation  (E.36)  is  correct.  Combining 
equations  (E.35),  (E.36),  (E.37),  (E.38),  and  (E.39),  and  using  the  fact 
that  fr  =  L/D,  we  obtain  finally  equations  (4.9)  and  (4.10)  in  Chapter 
4. 


Constraint  Boundaries  for  s-j  and  klfn.  From  physical  considerations 

we  require  a  nonpositive  slcpe  at  X-j  and  nonpositive  curvature  Kx  at 

the  maximum  diameter  section.  We  have  defined  Si  and  k,  so  that 

1  *m 

si  >  0  and  klm  >.  0  satisfy  these  physical  requirements. 

We  also  require  that  no  inflection  points  occur  on  the  midbody 
except  at  Xi.  The  limiting  inflection  point  occurs  when  the  second 
derivative  touches  zero  but  does  not  change  sign.  This  is  expressed  by 
equations  (E.19)  which  reduce  immediately  to 

y"(x)  =  0  (E.40a) 

y"’(x)  =  0  (E  .40b) 

Equations  (E.4Q)  can  be  solved  by  direct  substitution  to  obtain 


t  _  G"(x)  FSXx)  -  G"’(x)  F 2 ( X ) 
l«n  Fi(x)  F^x)  -  Fi^x)  Fg'(x) 


(E.41a) 


klm  Fi'(x)  -  G"(x) 

F^Txl 


( E .41 b) 


for  the  limiting  inflection  point  condition.  The  solution  is  tabulated 
in  Table  E3  and  plotted  in  Figure  21  in  Chapter  4. 
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There  is  a  singularity  in  equations  (E.41)  for  x  =  0.  For  this 
case  equation  (E.40a)  is  identically  satisfied  and  equation  (E.40b) 
reduces  to 


This  straight  line  is  plotted  in  Figure  23  in  Chapter  4. 


Tail  boom  Aftbody  Section 


The  tailboom  aftbody  section  used  with  the  eight-parameter  tail- 
boom  body  is  shown  in  Figure  E4  with  its  dimensional  parameters.  An 
inflection  point  occurs  at  X-j;  at  this  point  the  slope  is  S-j  and  the 
radius  is  R^.  There  is  also  an  inflection  point  at  L;  the  terminal 
radius  is  T.  The  boundary  conditions  for  this  profile  Y(X)  are  listed 


Figure  E4.  Tailboom  Aftbody  Section. 


below  where  primes  on  Y(X)  denote  differentiation  with  respect  to  X: 


W 'A KW.  RIUW* L F.fWI w 1  m-  W'»":^l  W '/JJVjlMlMpj  ■*-»> *>S'," 


Table  E3.  Solution  for  Limiting  Inflection  on  Midbody. 


X 

V  * 

-m 

si 

X 

k  * 

% 

si 

.020 

5.0511 

1.2457 

.520 

7.6460 

1.2787 

.04G 

5.1048 

1.2415 

.540 

7.7794 

1.3161 

.050 

5.1612 

1.2371 

.560 

7.8913 

1 .3641 

.080 

5.2203 

1.2327 

.580 

7.9698 

1 .4246 

.100 

5.2826 

1.2282 

.600 

8.0000 

1.5000 

.120 

5.3480 

1  .2237 

.620 

7.9632 

1.5919 

.140 

5.4169 

1.2192 

.640 

7.8387 

1.7016 

,160 

5.4896 

1.2147 

.660 

7.6052 

1 .8289 

.180 

5.5662 

1.2102 

.680 

7.2452 

1.9716 

.200 

5.6470 

1.2058 

.700 

6.7500 

2.1250 

.220 

5.7324 

1.2016 

.720 

6.1250 

2.2812 

.240 

5.8225 

1.1975 

.740 

5.3936 

2  4308 

.260 

5.9178 

1.1938 

.760 

4.5957 

2.5638 

.280 

6.0185 

1.1904 

,780 

3.7812 

2.6718 

.300 

6.1250 

1.1875 

.800 

3.0000 

2.7500 

.320 

6. 7374 

1.1852 

.820 

2.2924 

2.7971 

.340 

6.356G 

1.1838 

.840 

1.6842 

2.8157 

.360 

6 . 48 1 0 

1.1835 

.860 

1.1854 

2.8104 

.380 

6.6123 

1.1846 

880 

.7941 

2.7357 

.400 

6.7500 

1.1875 

.900 

.5000 

2.7500 

.420 

6.8934 

1.1926 

.920 

.2891 

2.7048 

.440 

7.0419 

1.2005 

.940 

.1467 

2.6548 

.460 

7,1940 

1.2121 

.960 

.0588 

2.6029 

.480 

> .3478 

1.2282 

.900 

.0132 

2.5508 

.500 

7.5000 

1.2500 

1.000 

.0000 

2.5000 

*  kim  =  (xi/xm  -  1 )2/ (1  ~  ri)  by  equation  (E.35) 
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1.  Y(X  j )  =  R-j  (E.43a) 

2.  Y 1 ( X i )  =  Si  (E.43b) 

3.  r(Xi)  =  0  (E.43c) 

4.  Y(L)  =  T  ( E . 43d ) 

5.  Y'(L)  =  0  (E,43e) 

6.  Y"(L)  =  0  (E.^oif ) 


where  Rj,  Sj,  T,  X-j ,  and  L  may  all  vary. 

It  is  possible  to  reduce  Y(X)  to  a  nondimensional  profile  with  only 
two  variable  boundary  conditions.  We  define  the  nondimensional  profile 


to  be 

x  =  (L  -  X)/(L  -  Xi)  (E.44a) 

y  =  Y/R-j  (E  .44b; 

so  that  the  boundary  conditions  given  by  equations  (E.31)  become 

1.  y(0)  =  T/R-,-  =  t/n  (E.45a) 

2.  y'(0)  =  0  (E .45b) 

3.  y"{0)  =  0  (E.45c) 

4.  y(l)  =  1  (E.45d) 

5.  y'  (1 )  =  -(L  -  X,-)  Si/Ri  =  Si  (E.45e) 

d 


6.  y"(l ) 


0 


(E.45f) 
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where  primes  on  y(x)  denote  differentiation  with  respect  to  x  and  t  = 
2T/D  and  r-j  =  2R-j/D.  It  is  apparent  from  equations  ( E . 34b)  and  (E.45e) 
that 

sia  =  s-j [ ( L  -  X-j ) (D/2  -  Ri)]/[(Xi  -  Xm) (R-j )3 

=  Si[(l  -  xi ; (1  -  ri)]/[(xi  -  xm)(r-j)]  (E .46) 

where  xm  =  Xm/L,  x-j  =  X-j/L,  and  r^  =  2R1-/D.  The  variable  boundary  con¬ 
ditions  are  sn-  and  t/r-,-;  two  free  parameter:  make  it  possible  to  con- 

1  a  1 

veniently  plot  constraint  boundaries  for  the  tailboom  aftbody. 

Following  the  procedure  in  Section  E.3,  we  postulate  immediately 

that 

(x)  -  ^  Fj(x)  +  Sja  F2(x)  +  G(x)  (E.47) 

where  Fx(x),  F2(x),  and  G(x)  are  polynomials  of  degree  five  since  there 
are  six  boundary  conditions  given  by  equations  (E.45).  For  arbitrary 
values  of  t/r-j  and  s-ja  the  boundary  conditions  for  the  function  Fj(x), 

F2(x),  and  G(x)  are 

F,(0)  =  1  F^l)  =  0 

FJ(0)  =  0  Fi(l)  =  0 

F”(0)  =  0  FJd)  =  0 

F2(0)  =  0  F2(l)  =  0 

F 1(0)  =  0  F^(l)  =  1 

f;(o)  =  o  F^(i)  =  o 

G(0)  =  0  G(1 )  =  1 

G'(0)  =  0  G'(l)  = 

G"{0)  =  0  G"d )  -  J 
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For  the  function  Fj(x)  we  may  write 

fa  Fx(x)  =  CjX2(x  -  l)2 

Integrating  and  applying  conditions  F^O)  =  1  and  Fx (1 )  =  0  gives 

Fj (x)  =  1  -  x3(6x2  -  15x  +  10)  (E.48) 

For  the  function  F2(x)  we  may  write 

5 

F2(x)  -  l  bj  x1 

i=0 

since  the  boundary  conditions  cannot  be  satisfied  by  the  simple  form 
(constant)  xm  (x  -  l)n.  By  direct  elimination  we  obtain 

F2(x)  =  -::3(3x2  -  7x  +  4)  (E.45) 

The  function  6(x)  has  the  same  bcndary  conditions  as  in  Section  E.3; 
thus  equation  (E.39)  applies  directly.  Equations  (E.39)  and  (E.48) 
imply  that 

G(x)  =  1  -  Fj (x)  (E.50) 

so  that 

y(x)  =  1  '  (77-  F,(x)  +  5ja  F2(x)  (E.51) 

Equations  (E.48),  (E.49),  and  (E.51)  satisfy  the  boundary  conditions 
(E.45);  hence  the  postulated  form  given  by  (£.47)  is  correct.  The 
results  are  given  as  equations  (4.11)  and  (4.12)  in  Chapter  4,  to  which 
equation  (E.46)  has  been  applied. 
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Constraint  Boundaries  for  s.;a  and  t/r-j .  From  physical  considera¬ 
tions  we  require  a  nonpositive  slope  S-j,  a  positive  profile  radius  Rj 
at  X-j,  and  a  non-negative  profile  radius  T  at  We  have  defined  s-ja> 

t,  and  r^  so  that  s-ja  0  and  0  <.  t/r^  <  1  satisfy  these  physical 
requirements. 

We  also  require  that  no  inflection  points  occur  on  th-4  tailboom 
aftbody  other  than  at  X-j  and  L;  this  implies  that  at  most  one  inflec¬ 
tion  can  occur  on  X -j  <  X  <  L.  Equations  (E .40)  cannot  be  used  to  find 
the  constraint  boundaries;  no  solution  exists  on  0  <  x  <  1  so  that  an 
alternate  approach  is  required.  Using  equation  (E.40a)  together  with 
equation  ( E . 51 )  one  may  write 

y"(x)  =  (£•-  1)  F;(x)  +  Sia  F^(x)  =  0  (E.52) 


Factoring  out  the  solutions  at  x  ^  and  x  =  1 ,  we  obtain  after  some 
manipulation  that 


x 


which  is  the  location  of  the  third  inflection  point.  It  is  apparent 
that  we  must  find  the  range  of  s-}a/ ( t/r-j  -  1)  for  whi  ,h  x  <  0  and  x  >.  1 . 

Setting  x  <  0  and  then  x  _>  1  in  equation  (E.53)  and  performing  the  in¬ 
equality  analyses  gives 


(E.54) 


which  defines  the  region  for  which  no  inflection  point  exists  on  0  <  x  < 
1.  These  linear  boundaries  are  plotted  in  Figure  24  of  Chapter  4. 
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APPENDIX  F 
ADDITIONAL  RESULTS 


This  appendix  contains  additional  optimization  results  not  reported 
in  Chapter  5.  Also,  summary  tablts  and  figures  of  all  optimization  runs 
are  presented.  Unless  otherwise  noted,  all  optimization  runs  are  made 
using  30-station  solutions.  However,  the  reported  drag  coefficients  are 
99-station  solutions.  The  computer  running  time  in  minutes  is  nominally 
equal  to  the  number  of  performance  function  (PF)  evaluations  on  the 
CDC  6500. 

F.l  Eight-Parameter  Tailboom  Body  at  Ry  =  7  x  106. 


The  results  for  the  eight-parameter  tailboom  body  are  limited;  the 
summary  is  shown  in  Table  FI.  The  additional  data  not  previously  re¬ 
ported  is  the  comparison  of  the  Cp  values  based  on  the  30-station  and 
99-station  solutions.  For  body  D-54  at  Ry  -  7  x  106  the  predicted  Cp 
decreases  112  when  99  stations  are  used. 


F.2  Five-Parameter  All-turbulent  Body  at  Ry  =  5  x  106. 


The  three  optimization  runs  for  this  body  are  summarized  in  Table 
F2;  overall  parameter  migrations  are  shown  in  Figure  FI.  From  the  table 
it  is  evident  that  there  is  essentially  no  relative  distortion  of  the 
response  surface  due  to  the  30-station  solutions.  This  suggests  that 
the  coarse  30-station  grid  loses  important  information  associated  with 
the  rapid  change  in  the  boundary  layer  immediately  downstream  of  the 
point  transition.  This  rapid  change  does  not  occur  in  an  all-turbulent 
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0—0  Series  !  (.01976) 

D— O  Series  12  (.02065) 

A— A  series  13  (.01982) 

Best  C~  of  Series  Shown  in  Parentheses 
i  0 

i 


i 


Figure  FI.  Summary  of  Parameter  Migrations  for 
All -turbulent  Bodies  at  Ry  =  5  x  10*. 
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boundary  layer.  From  Figure  FI  it  is  apparent  that  the  response  surface 
for  all-turbulent  bodies  at  Ry  =  5  x  106  is  quite  flat.  The  three  solu¬ 
tions  are  spatially  far  apart  but  their  Cq  values  are  all  within  4.5% 
of  the  best  value.  The  parameters  for  the  three  all-turbulent  solutions 
are  given  in  Table  F3. 

F.3  Five-Parameter  Laminar  Body  at  Ry  =  5  x  106. 

The  optimization  runs  for  this  body  are  summarized  in  Table  F4; 
overall  parameter  migrations  are  given  in  Figure  F2.  It  is  evident  in 
the  table  that  there  is  significant  relative  distortion  of  the  response 
surface  due  to  the  30-station  solutions.  The  Complex  Method,  moving  on 
global  information,  can  cope  with  the  relative  distortion;  that  is,  the 
distortion  apparently  does  not  obliterate  global  trends.  The  Powell 
Method,  moving  on  local  information,  is  useless  when  using  30-station 
solutions;  the  results  of  optimization  Series  P  in  Table  F4  demonstrate 
this.  The  results  of  Series  P  motivated  the  initiation  of  Series  T 
which  uses  99-station  sclutir.is.  Of  course,  this  procedure  triples 
the  cost  of  the  op  4  rotation  run. 

F.  t  Five  Parameter  Laminar  Body  at  Ry  =  1.6  x  IQ7. 

The  optimization  runs  for  this  body  are  summarized  in  Table  F5; 
overall  parameter  migrations  arc  shown  in  Figure  F3.  As  discussed 
above,  there  is  relative  distortion  of  the  response  surface  due  to  the 
30-station  solutions.  It  appears  that  99-station  solutions  must  be 
cm pi oyed  to  use  the  Powell  Method  effectively. 


Table  F3.  Parameter  Values  for  Three  All -turbulent  Bodies 
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F.5  Five-Parameter  Laminar  Body  at  P.y  =  5  x  107. 


Extensive  runs  have  been  made  at  this  Reynolds  number  in  an  attempt 
to  establish  the  existence  of  uniqueness  of  the  minimum  drag  shapes. 

These  runs  are  summarized  in  Table  F6;  the  overall  parameter  migrations 
are  shown  in  Figure  F4.  The  results  show  that  response  surface  distor¬ 
tion  is  present  as  mentioned  in  the  preceding  cases;  however,  it  appears 
that  large  scale  trends  are  not  obliterated. 

Based  on  f  *ies  L,  far  from  the  feasible  optimum  the  response  sur¬ 
face  is  nominally  flat  with  local  minima  present.  The  series  is  termi¬ 
nated  before  formal  convergence  is  achieved,  but  after  70  PF  evaluations 
the  method  is  makinq  little  progress.  The  three  Powell  series  K,  M,  and 
N  all  move  rapidly  co  the  rn  -  kx  boundary  and  then  move  slowly  along 
it.  It  appears  that  the  Powell  Method  used  in  conjunction  with  the 
penalty  function  defined  by  equation  (3.13)  cannot  cope  effectively  with 
curved  boundaries.  Stated  another  way,  the  direction-seeking  strategy 
of  the  method  has  difficulty  aligning  the  search  along  boundaries  which 
are  not  parallel  t.o  the  initial  set  of  search  directions,  which  usually 
is  the  set  of  parameter  axes. 

Series  F  and  F2,  both  using  the  Complex  Method,  converge  to  distinct¬ 
ly  different  shapes,  as  reported  in  Chapter  5.  For  this  reason,  the 
region  containing  these  two  solutions  is  examined  in  detail  to  determine 
whether  the  solutions  are  unique.  By  uniqueness  we  mean  the  existence 
of  a  finite  number  of  distinct  minima.  Series  R,  using  the  Complex 
Method,  revealed  no  significant  information  regarding  uniqueness. 

The  best  body  shape  obtained  formally  by  an  optimization  run  is 
body  M-73  which  lies  on  the  rn  -  kj  boundary.  Mainly  due  to  the  trend 
of  Series  F2  and  M,  it  is  believed  that  the  minimum  Cg  at  this  Reynolds 
number  lies  on  this  boundary.  Perturbations  on  body  M-73  along  the 
rn  -  kx  boundary  are  shown  in  Figure  F5.  A  generally  well  behaved  trend 
is  evident;  the  bracketed  minimum  is  near  solution  MP7  with  a  Co  value 
of  .00712.  This  behavior  suggests  that  the  minimum  drag  solution  is 
unique.  However,  this  cannot  be  construed  as  a  sufficient  test  for 
uniqueness;  a  truly  sufficient  test  for  uniqueness  does  not  exist  for 
this  problem. 
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Figure  F5.  Perturbations  on  Body  M-73  along  rn  -  Boundary. 


Assuming  that  the  optimum  solution  at  Ry  =  5  x  107  is  body  MP7,  then 
it  is  proper  to  question  the  validity  of  the  presumed  boundary  on  which 
the  optimum  lies.  An  additional  perturbation,  body  MP5  with  Cq  *  .01279, 
reveals  that  the  response  surface  ascends  rapidly  outside  the  presumed 
feasible  region.  This  tends  to  confirm  the  previous  hydrodynamic  expe¬ 
rience  which  motivated  the  presence  of  this  boundary  at  the  outset.  The 
parameters  for  six  laminar  bodies  obtained  in  studies  at  Ry  =  5  x  107 
are  summarized  in  ’’"able  F7. 

It  is  of  interest  to  examine  the  hydrodynamic  reasons  why  body  MP7 
has  a  low  Cq  while  its  close  neighbor,  body  MP5,  has  such  an  inferior 
Cq  value.  Body  MP7  and  its  velocity  distribution  are  shown  in  Figure 
F6.  The  ^orebody  has  a  limiting  inflection  at  X/L  *  .26.  The  curvature 
behavior  induces  an  early  adverse  velocity  gradient  region  which  the 
laminar  boundary  survives  at  this  Reynolds  number,  at  least  according  to 
the  flow  model.  A  region  of  locally  accelerated  flow  suppresses  transi¬ 
tion  until  X/L  =  0.5.  The  early  adverse  region  helps  to  reduce  skin 
friction  but  it  has  a  destabilizing  effect  on  the  laminar  flow.  This 
can  be  seen  in  Figure  F7  which  shows  the  Rq  versus  R$  trajectory  at 
Ry  =  5  x  107.  The  trajectory  very  nearly  touches  the  Michel -e’  curve 
at  X/L  *  .25  approximately.  The  region  of  locally  accelerated  flow 
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Table  F7.  Parameter  Values  for  Six  Laminar 
Bodies  Designed  at  Ry  *  5  x  107. 


Parameter 

F-57 

F2-49 

Body  Number 

M-73  MP5 

MP7 

R-2 

*r 

4.2735 

3.5000 

3.5306 

3.5306 

3.5306 

3.7866 

*m 

.4446 

.4300 

.4710 

.4710 

.4710 

.4729 

3.8081 

3.7000 

3.8191 

3.5000 

3.3500 

4.1106 

?r\ 

.1821 

.9000 

1.0217 

1.5500 

1 .3900 

.5607 

A 

.1773 

.2000 

.1789 

.1789 

.1789 

.0478 

*  k>a 

5.9443 

6.5015 

4.8157 

4.4133 

4.2242 

5.1086 

CD 

.0077 

.0073 

.0072 

.0128 

.0071 

.0072 

ty  Ratio 


183 


Log  R$  =  Log  S  ue/v 
Figure  F7.  Re  Versus  R$  for  Body  MP7. 


causes  the  trajectory  to  veer  away  and  then  cross  tie  correlation  curve 
at  X/L  *  0.3.  The  neighboring  solution,  body  MP5,  appears  very  similar 
to  body  MP7  and  has  a  similar  velocity  distribution  (figure  not  shown). 
However,  the  early  adverse  region  for  body  MP5  is  severe  enough  to 
cause  its  R@  versus  R5  trajectory  to  cross  the  Michel-e*  curve  at  X/L  * 
.22.  This  accounts  for  the  rapid  change  in  Cp  between  the  neighboring 
solutions. 

The  high  sensitivity  of  body  MP7  to  early  transition,  as  Inferred 
from  Figure  F7,  makes  it  an  undesirable  hydrodynamic  design.  The  trade¬ 
off  between  low  Cq  and  low  sensitivity  to  early  transition  must  be  left 
to  the  judgement  of  the  designer.  It  is  possible,  of  course,  to  Insert 
additional  constraints  into  the  drag  minimization  problem  to  avoid  this 
undesirable  sensitivity. 


